\documentclass[12pt,Bold,letterpaper]{mcgilletdclass}
\usepackage[dvips,final]{graphicx}
\usepackage[top=1in,bottom=1in,dvips]{geometry}
\usepackage{float}		% this is to place figures where requested!
\usepackage{times}		% this uses fonts which will look nice in PDF
\usepackage{graphicx}		% needed for the figures
\usepackage{epstopdf}
\usepackage{url}
\usepackage{amsfonts }
\usepackage{amssymb,amsmath}
\usepackage{algorithm}
\usepackage{verbatim}
\usepackage{algorithmic}
%\usepackage{program}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{relsize}
\usepackage[english,french]{babel}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\newcommand{\bmx}{\begin{bmatrix}}
\newcommand{\emx}{\end{bmatrix}}
\newcommand{\bsmx}{\left[\begin{smallmatrix}}
\newcommand{\esmx}{\end{smallmatrix}\right]}

% Equations:
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\beqy}{\begin{eqnarray}}
\newcommand{\eeqy}{\end{eqnarray}}
\newcommand{\beqynn}{\begin{eqnarray*}}
\newcommand{\eeqynn}{\end{eqnarray*}}
\newcommand{\boxcon}{{\cal B}}
\newcommand{\bd}{{\bar{d}}}
\newcommand{\bl}{{\bar{l}}}
\newcommand{\bu}{{\bar{u}}}
\newcommand{\br}{{\bar{r}}}
\newcommand{\bx}{{\bar{x}}}
\newcommand{\by}{{\bar{y}}}
\newcommand{\bz}{{\bar{z}}}
\newcommand{\dist}{\mathrm{dist}}
\newcommand{\vsp}{\vspace{\baselineskip}}

%\usepackage{tex4ht}
%\usepackage{amsmath}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Have you configured your TeX system for proper  %%
%% page alignment? See the McGillETD documentation %%
%% for two methods that can be used to control     %%
%% page alignment. One method is demonstrated      %%
%% below. See documentation and the ufalign.tex    %%
%% file for instructions on how to adjust these    %%
%% parameters.                                     %%
\addtolength{\hoffset}{0pt}                        %%
\addtolength{\voffset}{0pt}                       %%
%%                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       Define student-specific info
\SetTitle{\huge{Integer Least Squares\\Search and Reduction Strategies}}%
\SetAuthor{Stephen Breen}%
\SetDegreeType{Master of Science}%
\SetDepartment{School of Computer Science}%
\SetUniversity{McGill University}%
\SetUniversityAddr{Montreal,Quebec}%
\SetThesisDate{August, 2011}%
\SetRequirements{A thesis submitted to McGill University\\
in partial fulfilment of the requirements of the degree of\\
Master of Science in Computer Science}%
\SetCopyright{\copyright Stephen Breen 2011}%

\makeindex[keylist]
\makeindex[abbr]

%% Input any special commands below
%\newcommand{\Kron}[1]{\ensuremath{\delta_{K}\left(#1\right)}}
\listfiles%
\begin{document}

\maketitle%
\justifying
\begin{romanPagenumber}{2}%

\SetDedicationName{\MakeUppercase{Dedication}}%
\SetDedicationText{\em \center To my patient and supportive parents Gwendolyn and Stephen Breen}%
\Dedication%

\SetAcknowledgeName{\MakeUppercase{Acknowledgements}}%
\SetAcknowledgeText{
First I would like to thank my supervisor, Xiao-Wen Chang for his excellent guidance and many helpful comments and ideas. From him I have learned the meaning of being rigorous and precise in ones work and how imprecise I often can be.\\
\indent I would also like to thank NSERC and McGill University for their generous financial support, without which this thesis would not have been possible.\\
\indent Thanks to my lab mates in the Scientific Computing lab, especially Wen-Yang Ku, Xiaohu Xie and Sevan Hanssian for the numerous times that they have helped me through difficult problems.\\
\indent I would like to thank my co-workers, Todd Murphy, Kevin Garnett, Jackie Howatt, and Allan Cameron at MedRunner, who have been very understanding with my lack of availability over the past two years.\\
\indent Finally I thank my friends and family who have remained close over time and distance. Especially my good friend Marya Sawaf to whom I owe the excellent last minute editing and French translation of the abstract. Thanks to her for always being there.
}%
\Acknowledge%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         English Abstract                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\SetAbstractEnName{\MakeUppercase{Abstract}}%
\SetAbstractEnText{
This thesis is concerned with integer least squares problems, also referred to as closest vector problems. One often used approach to solving these problems is the discrete search method, which typically involves two stages, the reduction and the search. The main purpose of the reduction is to make the search faster.\\
\indent Reduction strategies for box-constrained integer least squares problems involve column reordering of the input matrix. There are currently two algorithms for column reordering that are most effective for the search stage, referred to here as SW and CH. Although both use all available information in the problem, the SW and CH algorithms look different and were derived respectively from geometric and algebraic points of view. In this thesis we modify the SW algorithm to make it more computationally efficient and easier to comprehend. We then prove that the SW and CH algorithms actually give the same column reordering in theory. Finally, we propose a new mathematically equivalent algorithm, which is more computationally efficient and is still easy to understand. This thesis also extends the column permutation idea to ordinary integer least squares problems. A new reduction algorithm which combines the well-known Lenstra–Lenstra–Lovász (LLL) reduction and the new column reordering strategy is proposed. The new reduction can be much more effective than the LLL reduction in some cases.\\
\indent The thesis also reviews some common search algorithms. A new one is proposed, which is based on two previous algorithms, the depth-first search and the best-first search. This hybrid algorithm makes use of the advantages of both originals, is more efficient than either and is easier to implement than other previous hybrid algorithms.
}
\AbstractEn%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         French Abstract                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\SetAbstractFrName{\MakeUppercase{ABR\'{E}G\'{E}}}%
\SetAbstractFrText{Cette thèse s'intéresse aux problèmes de moindres carrés entiers (ILS), ou les problèmes du vecteur le plus proche. Une approche souvent utilisée pour résoudre ces problèmes est la méthode de recherche discrète, qui implique deux étapes: la réduction et la recherche. Le but principal de la réduction est de rendre l’étape de recherche plus rapide.\\
\indent Les stratégies de réduction des problèmes ILS sous contrainte de boîte impliquent la réorganisation de colonnes de la matrice de données. Il existe actuellement deux algorithmes pour la réorganisation des colonnes, appelés ici les algorithmes SW et CH,  qui sont les plus efficaces pour la phase de recherche. Bien que les deux utilisent toutes les informations disponibles dans le problème, les algorithmes SW et CH sont différents en apparence, et ont été obtenus respectivement à partir d’une point de vue géométrique et algébrique de vue. Dans cette thèse, nous modifions l'algorithme SW pour rendre son calcul plus efficace et plus facile à comprendre. Nous démontrons ensuite qu’en théorie, les algorithmes SW et CH donne effectivement la même réorganisation de colonnes. Enfin, nous proposons un nouveau algorithme mathématiquement équivalent qui est plus efficace, tout en demeurant facile à comprendre. Cette thèse étend également l'idée de permutation de colonnes aux problèmes ordinaires de moindres carrés entiers. Un nouveau algorithme de réduction qui combine le célèbre agorithme Lenstra-Lenstra-Lovász (LLL) avec la nouvelle stratégie de réorganisation de colonnes est proposé. La nouvelle réduction peut être plus efficace que la réduction LLL dans certains cas.\\
\indent Cette thèse examine également certains algorithmes de recherche d’usage courant. Un nouveau est proposé qui est basé sur deux algorithmes précédents: l’algorithme de parcours en profondeur et celui de la recherche au meilleur d'abord. Notre algorithme hybride détient les avantages des deux originaux, tout en étant plus efficace et plus facile à utiliser que d’autres algorithmes hybrides déjà existants.
}%
\AbstractFr%

\TOCHeading{\MakeUppercase{Table of Contents}}%
\LOTHeading{\MakeUppercase{List of Tables}}%
\LOFHeading{\MakeUppercase{List of Figures}}%
\tableofcontents %
\listoftables %
\listoffigures %

\end{romanPagenumber}

%\mainmatter %
\selectlanguage{english}

\chapter{Introduction}

This thesis deals with efficient algorithms for solving integer least squares (ILS) problems, sometimes referred to as closest vector problems (CVP). This problem has been well studied in mathematics for over $100$ years \cite{HanPS11}. The ordinary ILS problem has been proven to be NP hard \cite{Mic01}, meaning that the time needed to solve the problem can depend exponentially on the problems size. Since ILS problems arise in a number of time sensitive practical applications, it is worthwhile to study algorithms which can solve them quickly. In this thesis an efficient class of algorithms for solving ILS problems will be studied. New algorithms will also be proposed which offer decreased computational or memory cost when used to solve problems that arise in some practical applications.

In this chapter an introduction to ordinary real least squares and integer least squares problems will be given; the distinction between the two will be made clear. Some applications for the ILS problem will be mentioned and then previous work which has dealt with solving these problems efficiently will be reviewed. 

\section{Notation}

This section introduces the notation which will be used throughout this thesis. Matrices will be denoted by capital letters while vectors and scalars will be in lower case. Subscripts will be used to identify elements of vectors and matrices and also columns of matrices. For example, if we have a matrix $A$, $a_j$ denotes the $j^{th}$ column and $a_{i,j}$ denotes the element in row $i$, column $j$. For a vector $v$, $v_j$ denotes the $j^{th}$ element. The identity matrix of size $n$ will be denoted as $I_n$ or $I$ and the $i^{th}$ column of the identity matrix is denoted as the vector $e_i$. Occasionally MatLab like notation will be used to denote subsets of rows/columns in a matrix or vector. For example $A_{:,j:k}$ denotes the matrix formed by columns $j, j+1, \dots, k$ of the matrix $A$.

Rounding of a real valued scalar or the elements of a real valued vector to the nearest integer(s) will be written as $\lfloor x \rceil$, where $x$ is the scalar or vector. Suppose we have some set of integers ${\cal B}$, rounding of $x$ to the nearest value in this set will be denoted as  $\lfloor x \rceil_{{\cal B}}$.

For a random vector $v$ following a normal distribution with $0$ mean and covariance matrix $\sigma^2I$, we write $v \sim N(0,\sigma^2I)$.

\section{List of Abbreviations}

The following is a list of abbreviations which will often be used in this thesis:
\begin{itemize}
\item LS: Least Squares.
\item ILS: Integer Least Squares.
\item RLS: Real Least Squares.
\item OILS: Ordinary Integer Least Squares (no constraints on the solution).
\item BILS: Box constrained Integer Least Squares (box constraint on the solution).
\item MIMO: Multiple Input, Multiple Output.
\item SE: The Schnor-Euchner Enumeration search strategy for ILS problems, \cite{SchE94}.
\item BFS: Best First Search, a tree search algorithm.
\item DFS: Depth First Search, another tree search algorithm.
\item CH: Chang and Han's algorithm for BILS reduction.
\item SW: Su and Wassell's algorithm for BILS reduction.
\item LLL: The Lenstra-Lenstra-Lovasz reduction strategy.
\item IGT: Integer Gauss Transformation
\end{itemize}

\section{Ordinary Real Least Squares Problem}
Consider the following linear model,
\begin{equation}
\label{eq:realLSModel}
y = Hx+v,
\end{equation}
where $y\in\mathbb{R}^m$ is the observation vector, $H\in\mathbb{R}^{m \times n}$ is called the
``design matrix'' and has full column rank, $x \in \mathbb{R}^n$ is the unknown parameter vector and $v \in \mathbb{R}^m$ is the noise vector with $v \sim N(0,\sigma^2I)$. To estimate $x$ from \eqref{eq:realLSModel}, an often used approach is to solve the following ordinary real least squares problem:
\be
\label{eq:realLSResidual}
	\min_{x \in \mathbb{R}^n} \left \| y- Hx \right \|^2_2.
\ee

The solution to this problem, $x_{\scriptscriptstyle{\mathrm{RLS}}}$ has the property that it is the maximum likelihood estimator of the true parameter vector $x$. This means that $x_{\scriptscriptstyle{\mathrm{RLS}}}$ is more likely than any other vector to be equal to the true, unknown parameter vector $x$. 

If we expand equation \eqref{eq:realLSResidual} and set its gradient to $0$, we will arrive at the
well known ``normal equations'', which can be written in matrix form as,
\begin{equation}
H^THx = H^Ty.
\label{eq:normalEquations}
\end{equation}
The solution to these normal equations satisfies the following:
$$x_{\scriptscriptstyle{\mathrm{RLS}}} = (H^TH)^{-1}H^Ty.$$
Some numerical methods to compute the solution of these least squares problems efficiently and reliably can be found in \cite{Bjo96}.

Least squares estimation has numerous applications in many fields in science, engineering and economics.

\vsp \section{Integer Least Squares Problems}
In some applications, the unknown parameter vector $x$ in equation \eqref{eq:realLSModel} is an integer vector, to estimate $x$ in this case, we instead solve the following integer least squares (ILS) problem:

\be
\label{eq:ils0}
\min_{x \in {\mathbb{Z}}^n}  \| y- Hx \|_2. 
\ee

We no longer have a closed-form representation for the optimal solution $x_{\scriptscriptstyle{\mathrm{ILS}}}$ in this case, in fact, the
problem is provably NP-hard \cite{Mic01}.  

Sometimes the integer parameter vector $x$ is subject to constraints. One example of such constraints is the following set of box constraints:
\begin{align} \label{eq:boxCon}
&x\in \boxcon,\\
&{\cal B} = {\cal B}_1  \times \cdots \times {\cal B}_n,\\
&\boxcon_i = \{x_i\in\mathbb{Z} : l_i\le x_i \le u_i, l_i\in \mathbb{Z},
u_i\in\mathbb{Z} \}.
\end{align}
The constrained integer least squares problem corresponding to these constraints is referred to as the box-constrained integer least squares (BILS) problem. 

Even though the problem \eqref{eq:ils0} is NP-hard, we still have some hope to get solutions
quickly for certain instances. The problems that arise in many practical applications can often be solved in reasonable amounts of time if efficient algorithms are used.

In lattice theory, the matrix $H$ is called the ``lattice generator matrix''. This is because it generates a lattice ${\cal{L}}(H)$ defined as follows:
\begin{equation}
{\cal{L}}(H) = \{Hx: x\in \mathbb{Z}^n\}.
\end{equation}
The columns of $H$ are called the basis vectors of the lattice. A single lattice may have many different sets of basis vectors or equivalently, many different generator matrices $H$. Figure \ref{fig:latticeBasis} shows an example of a 2 dimensional lattice with two different sets of basis vectors.

\begin{figure}
\centering
\includegraphics[scale=0.4]{latticebasis.png}
\caption{An example of a lattice with two different sets of basis vectors.}
\label{fig:latticeBasis}
\end{figure}

In this context, the OILS and BILS problems can be thought of as trying to locate the closest lattice point to some real point $y$. This is why the problem is sometimes referred to as the closest vector problem (CVP).

\vsp \section{Applications}
One application of the OILS problem arises in GPS, where carrier phase measurements are used. In GPS, there are two types of measurements that can be used to determine the position of a receiver, code phase and carrier phase. Code phase measurements can give accuracy to a few meters, while carrier phase are accurate to centimeters. To make use of the more accurate carrier phase measurements, we must know how many cycles the carrier wave has gone through between a satellite and a receiver. GPS can only directly know the fractional part of the number of cycles and the difference in the number of cycles between two successive measurements. The problem of estimating the integer part of the number of cycles is called integer ambiguity estimation and one of the usual approaches taken to solve it is to construct and solve an ILS problem. For details on how the problem is formed and the usual approaches taken to solve this problem, see \cite{Xu07}. 

Some important applications such as MIMO (multiple-input, multiple-output) wireless signal decoding depend on the
solution of the BILS problem. MIMO refers to the case where a wireless system has multiple
input antennas transmitting a signal which is received by multiple output
antennas. The purpose of a MIMO system is to maximize the throughput of the communications channel. Throughput is defined as ``the average rate of successful message delivery over a communication channel''. 

The signal received is our input vector $y$ from
equation \eqref{eq:realLSModel}, it has undergone some linear transformation by the known
``channel matrix'' $H$ (design matrix) and some noise has been introduced during
the transmission. Originally, we know that each element of $x$ came from some
finite set of symbols that we may want to transmit or receive (we model
this property with $\boxcon$). The purpose of such a system is to maximize
throughput, however, the overall throughput of the system depends on how
quickly and accurately we can solve the BILS problem. As the number of antennas increases, the BILS problem becomes the bottleneck on the throughput. In order to obtain estimates for $x$ quickly we may use suboptimal solutions to the BILS problem, but under the assumption that the noise has $0$ mean and is
normally distributed, the BILS solution is more likely than any other possible
solution to be the true integer parameter vector $x$ \cite{Han11}. For this reason, we say that a receiver which is decoding transmissions using an algorithm which solves the BILS problem exactly achieves ``optimal performance''; performance refers to the likelihood that the vector found by the decoder is equal to the transmitted vector $x$. For more information on this problem refer to \cite{Jan04}.

Other applications of OILS and BILS include cryptography, lattice design and number theory. For more details on some of these applications see \cite{HanPS11}.

\vsp \section{Previous Work} \label{sec:prevWork}
Due to the important applications of the OILS and BILS problems, much work has
gone into solving them efficiently. There are three main families of algorithms to solve these problems exactly. One family uses ``voronoi cells'', and will not be considered in this thesis because they are known to be relatively
inefficient \cite{HanPS11}. Another family of algorithms can be referred to as ``real
relaxation based branch and bound'' methods or RRBB methods. Using these algorithms to solve ILS
problems is a relatively new idea and no comprehensive comparisons have been
done to study their efficiency, although some preliminary tests indicate they
can be quite efficient in some circumstances and inefficient in others \cite{Ku11}. The last main family of algorithms
that solve ILS problems exactly are ``discrete search methods''. This thesis
will focus on studying these discrete search methods since they are often used in applications. There have also been some algorithms proposed that yield fast,
approximate solutions to the problems, some giving statistical bounds on the
likelihood of error. These algorithms can sometimes be referred to as
``randomized'' algorithms. This thesis however will only deal with finding the
optimal solution to the problem. For more details on the randomized algorithms
and Voronoi cell algorithms, see \cite{HanPS11} and references therein.

The discrete search based approaches which will be studied in this thesis try to find the optimal solution by enumerating vectors in the search space one element at a time until the optimal solution is found. Discrete search algorithms typically involve two stages, the reduction and the search. The first stage is the reduction, it performs manipulations on the matrix $H$ and vector $y$ from \eqref{eq:ils0} in order to optimize the problem for the search stage. There are a special set of manipulations which may be performed so that the solution is not modified. The search is performed after the reduction. In the search stage we enumerate integer points until we find the optimal one.

\vsp \subsection{Reduction Strategies} \label{subsec:Reductions}

As mentioned previously, the purpose of the OILS and BILS reduction is to transform the problem to an equivalent, but easier to solve one. This section will review some existing reduction strategies.

One very effective reduction strategy is known as the Korkine-Zolotarev (KZ) 
reduction \cite{KorZ73}. Unfortunately the KZ reduction is prohibitively expensive to compute,
it requires solving a number of so called ``shortest vector problems'' which are suspected to be NP-hard. The only time this reduction
may be practical is when a very large number of input vectors $y$ must be
decoded for the same matrix $H$ \cite{AgrEVZ02}. The standard reduction algorithm used in practice for reducing unconstrained OILS problems is the ``Lenstra-Lenstra-Lovasz'' or LLL
reduction \cite{LenLL82}. It is usually much more efficient to compute than the
KZ reduction, and, in fact, any KZ reduced matrix is also LLL reduced. For more details
on how the LLL reduction and KZ reduction are related, again see
\cite{AgrEVZ02}. The LLL reduction has been proven to converge in a finite
number of steps and to have average case polynomial complexity when the basis
vectors are normally distributed inside the unit sphere in $\mathbb{R}^m$ \cite{LinMH10}. There
are also results that bound the complexity based on the condition number of the
input matrix. For more details on the algorithms, complexity, properties and applications for LLL
reduction, see \cite{WubSJM11} and \cite{NguV10}.

In \cite{XieCB11}, it was found that many of the operations used in the original LLL algorithm are not always required as they do not affect the search process when the SE (Schnorr-Euchner) algorithm is used. The new reduction that results from applying only a subset of the operations is called the partial LLL reduction.

Unfortunately, applying LLL or partial LLL reduction will change a box constraint to a much more complicated constraint which will make the search process more difficult. For the BILS problem, there are other reduction strategies that focus only on permuting the columns of the matrix $H$, while the LLL reduction performs some other operations that are described in Chapter \ref{chap:Reduction}. We can separate reduction algorithms which simply try and find some optimal permutation of the columns of the matrix $H$ into two categories, those that only use the information contained in the matrix $H$, and algorithms that use all of the information in $H$, $y$ and the box constraint.

Two algorithms that only use the information in $H$ are the ``Vertical-Bell Laboratories Layered Space-Time'' or V-BLAST column reordering \cite{FosGVW99} and the ``Sorted QR Decomposition'' or SQRD column reordering \cite{WubBRKK01}. An examination of these two algorithms reveals very similar motivations behind each.

Algorithms that use all of the information in $H$, $y$ and the box-constraint are a fairly new development. The first was proposed by Su and Wassell in 2005 \cite{SuW05}, and the second by Chang and Han in 2008 \cite{ChaH05}. Numerical results show that these algorithms can offer great improvements over the previous reductions that use only the information contained in $H$.

\vsp \subsection{Search Strategies}
In Chapter \ref{chap:SESearch}, it will be shown that the discrete search process is equivalent to a tree search problem. The tree in question has a depth of $n$; each node on a path from root to leaf (where we define leaves as nodes at depth $n$) represents a valid element in the solution vector. Therefore the tree has exponential width. To see an example of a subset of the nodes in such a tree, refer to Figure \ref{fig:treeSearch} on page \pageref{fig:treeSearch}. The most widely used algorithm for the OILS search process, the Schnorr Euchner (SE) enumeration \cite{SchE94}, can be thought of as a fairly straight forward depth first search in such a tree. The SE search is based on previous work such as the Pohst \cite{Pho81} and Kannan \cite{Kan83} enumeration algorithms. Other tree search algorithms may also be used to solve the problem with varying degrees of efficiency.

In the literature, some modifications to the best first tree search strategy have also been proposed to solve this problem. A few such algorithms are given in \cite{MurGDC06}, \cite{XuWZW04}, \cite{FukMU04}, \cite{StuBF07} and \cite{DaiY08}. When doing a tree search, the disadvantage of the best first approach is that the memory requirements can be exponential in the worst case and there is a significant overhead to visit each tree node. This is compared to the depth first search where the memory requirements are linear in the problem size $n$ and there is very little computational cost to visit a node in the tree. The advantage of the best first approach is that it is guaranteed to visit the least number of nodes in the tree. Some of the papers listed above propose a pure best first search, while others try to make some sort of trade off to achieve lower memory usage.

There have been some attempts to compare different discrete search algorithms for the OILS search process. One such paper is \cite{MurGDC06}. The authors here devise a common framework (based on a tree search) that many search algorithms can be described within and from there they can do a comparison on the estimated computational complexity of each. Unfortunately through this theoretical comparison, we can only relate the number of nodes in the search tree that are visited by various algorithms, this does not fully consider the amount of time processing each node which is often a computational bottleneck.

There have been other suggestions to improve the discrete search based algorithms as well. In \cite{StoVH08}, it is proposed that by using lower bounds on the residual from the optimal solution, we can shrink the search space (equivalently, prune the search tree). A few such lower bounds are given for special cases of the OILS problem and one for the general OILS and BILS problem. Unfortunately, the computational complexity of computing some of these bounds can be prohibitive since it adds to the processing that must be performed at each node in the tree. Overall it seems that these lower bounds do not offer a decrease in computational complexity during the search process if the search is performed on a reduced problem.

Another method that attempts to shrink the search space is given in \cite{SchFL09}, in which the authors propose a simple stopping criteria for the search process that in theory should allow it to terminate earlier. Unfortunately, the bound derived here is not tight enough to be useful in practice and is rarely or never satisfied. To make the bound more likely to be satisfied, they propose to increase it by some amount. Unfortunately increasing the bound in this way means we may terminate with a suboptimal solution; this thesis will not consider this since it is only concerned with the exact solution of the OILS and BILS problems. 

\vsp \section{Objectives and Contribution}

This thesis will study both the search and reduction steps of the discrete search based ILS algorithms. For BILS reduction algorithms, we are limited to performing only column permutations on the matrix $H$; the reason for this will be explained in Chapter \ref{chap:Reduction}. There are a few algorithms which calculate how we should permute the columns of $H$, some of which were briefly introduced in section \ref{sec:prevWork}. Two more recent developments, \cite{ChaH05} and \cite{SuW05}, use all of the information in the matrix $H$, vector $y$ and the box-constraint to calculate the permutations. These algorithms have shown excellent results. In this thesis it will be proven that these two algorithms are theoretically equivalent even though they look much different. Knowing that they are equivalent, we can use the best ideas from both to create a new reduction strategy that is faster than either of the originals and is numerically reliable (the faster of the two original algorithms is not). Another advantage of these algorithms being equivalent is that since one had a geometric motivation and the other was derived algebraically, we now have both geometric and algebraic justification for why the column orderings given by these algorithms should help speed up the search process. Also, the algorithm in \cite{SuW05} was derived through a geometric motivation and as such is described in terms of geometry in the original paper. Their geometric interpretation is difficult to understand, perhaps due to page limit; this thesis will provide an algebraic explanation for the algorithm presented in \cite{SuW05} and offer some insight to the original. A shorter version of this work has been accepted for presentation at the IEEE GLOBECOM 2011 conference \cite{BreC11}.

The motivation for the permutation based reduction strategies is not specific to the BILS problem; these reduction algorithms can also reduce the run time for OILS problems as well. However, the very effective LLL reduction provides better results than using the above permutation strategy alone. One way to think about what each type of algorithm is doing is, the LLL reduction finds a new set of shorter and more orthogonal basis vectors, while the permutation based reductions are just finding an ordering for the basis vectors that performs well in the search process. Consider Figure \ref{fig:latticeBasis} and suppose the original matrix $H$ has columns $v_1$ and $v_2$, the LLL reduced matrix may have columns $u_1$ and $u_2$. The permutation based reductions are limited to simply changing the ordering so that $v_1 = v_2$ and $v_2 = v_1$. It is known that shorter more orthogonal basis vectors are preferable in the search process. By first performing LLL reduction to get a good set of basis vectors, and then applying a permutation based reduction to reorder them, we can sometimes greatly improve the performance of the search process. In this thesis, the strategy of first applying the LLL reduction, and then column permutations will be explored.

With many different algorithms for both the reduction and search process, it is not always clear how they relate. In \cite{MurGDC06} the authors propose a tree search framework to relate some of the various discrete search algorithms. This tree search framework makes the trade off between memory usage and the number of nodes visited in the tree search clear. Also we can see that as the memory usage increases, the computational complexity of visiting a node does as well. This thesis will consider combining two efficient search strategies (the SE search and best first search) in a new way which is not considered by this framework. A parameter will control the influence of each individual algorithm. Using actual run time simulation results from a comparison of this combined algorithm with various settings of its parameter, we can see the strengths of each of the different approaches. This allows us to tune the new combined algorithm so that it outperforms either of the originals in some practical applications. This also provides some insight into the runtime performance of each algorithm individually.

\vsp \section{Outline}
The rest of the thesis will be organized as follows.

In Chapter \ref{chap:SESearch}, the Schnorr-Euchner (SE) enumeration algorithm \cite{SchE94} will be presented in detail, much of the remainder of the thesis will use ideas and notation which comes from this algorithm. Also, since the reduction processes are trying to optimize the search process, it is critical to first understand the search process before considering the reduction.

In Chapter \ref{chap:Reduction}, an explanation will be given for why we need different reduction strategies for OILS and BILS problems. Then, strategies for reducing OILS and BILS problems will be presented separately. After introducing two effective algorithms for BILS reduction from the literature, we will give a new interpretation and some improvements for one of them. Using this new interpretation, we can prove that the two algorithms are mathematically equivalent. A new BILS reduction strategy will then be introduced which is based on the equivalence of the two algorithms. Finally, we will introduce a new reduction strategy for OILS problems.

In Chapter \ref{chap:Searches}, some other notable search algorithms and modifications to the basic SE enumeration will be given. Also, a new hybrid search algorithm is proposed which combines two of the original algorithms in order to take advantage of the positive features of each.

Finally, Chapter \ref{chap:Conclusion} will give a summary and highlight problems which need further investigation.

\chapter{Schnorr-Euchner Enumeration} \label{chap:SESearch}

As mentioned previously, there are two steps to solve ILS problems by the discrete search approach, reduction and search. The purpose of the reduction step is to transform the problem to an easier but equivalent one so the search proceeds more quickly. Therefore the reduction is performed before the search. In this thesis the details on the search process are presented first for easier understanding. The reason for this is that the reduction manipulates the matrix $H$ in such a way to optimize it for the search process, therefore understanding the search process is crucial to understanding the reduction.

There are many search algorithms that have been proposed to solve the OILS problem. One of the most effective algorithms in terms of both overall runtime and memory consumption is the Schnorr-Euchner enumeration \cite{SchE94}. The SE search strategy is an improvement to a prior strategy called the Pohst enumeration strategy. For more details on the Pohst algorithm and how it relates to the SE algorithm, see \cite{AgrEVZ02}. In this chapter, the SE algorithm is presented in detail, since concepts from it will be used throughout the remainder of the thesis.

Let $H$ from the ILS problem in equation \eqref{eq:ils0} have the QR decomposition
$$
H=[Q_1, Q_2] \bmx R \\ 0 \emx,
$$
where $[\underset{n}{Q_1}, \underset{m-n}{Q_2}]  \in \mathbb{R}^{m\times m}$ is orthogonal
and $R\in \mathbb{R}^{n\times n}$ is upper triangular. 
Then, with $\bar{y}=Q_1^Ty$ the OILS problem \eqref{eq:ils0} is reduced to 
\be 
\label{eq:ils}
\min_{x \in  {\mathbb{Z}^n}}  \| \by- Rx \|_2.
\ee

The goal of the search strategy is to enumerate the values for the elements of $x$ until the optimal solution $x_{\scriptscriptstyle{\mathrm{ILS}}}$ is found. We would like to enumerate as few integer points as possible to find the optimal solution. Suppose that the optimal solution satisfies
\be 
\label{eq:searchIneq0}
 \| \by- Rx \|_2 < \beta. 
\ee
We will discuss some good ideas for choosing an initial $\beta$ later. The inequality \eqref{eq:searchIneq0} defines a hyper-ellipsoid in terms of $x$ or a hyper-sphere in terms of the lattice points $w=Rx$ with radius $\beta$. A search algorithm finds the optimal lattice point within this hyper-sphere. For this reason, the problem is sometimes referred to as ``Sphere Decoding'' in communications \cite{HasV05}. For an example of what the geometry may look like, refer to Figure \ref{fig:latticeBasis} and consider the two circles.

Define
\begin{equation}
 c_k = (\bar{y}_k - \sum_{j=k+1}^nr_{kj}x_j)/r_{kk}, \; k=n, n-1,\ldots, 1,
\label{eq:searchC}
\end{equation}
where when $k=n$ the sum in the right hand side does not exist. Notice that $c_k$ depends on $x_{k+1:n}$.
Then \eqref{eq:searchIneq0} can be rewritten as
\begin{equation}\label{eq:searchIneq1}
\sum_{k=1}^n r_{kk}^2(x_k-c_k)^2 < \beta^2,
\end{equation}
which implies the following
set of inequalities:
\begin{align}
&\text{level } k: \ \ r_{kk}^2(x_k-c_k)^2 < \beta^2 -\sum_{i=k+1}^nr_{ii}^2(x_i-c_i)^2, \label{eq:searchLevelK}
\end{align}
for $k=n,n-1,\ldots, 1$.

We begin the search process
at level $n$, therefore the valid values for $x_n$ only depend on the initial $\beta$. Choose $x_n = \lfloor c_n \rceil$, the nearest integer to $c_n$. If the inequality \eqref{eq:searchLevelK} with $k=n$
is not satisfied then it will not be satisfied for any integer. This means $\beta$
was chosen to be too small initially so it must be enlarged and the search restarted. With $x_n$ fixed, we can move
to level $n-1$ and choose $x_{n-1} = \lfloor c_{n-1} \rceil$ with $c_{n-1}$ calculated as in equation \eqref{eq:searchC}. At this point it is possible that the inequality \eqref{eq:searchLevelK} is no longer satisfied. If this is the case, we must move back to level $n$ and choose $x_n$ to be the second nearest integer to $c_n$.  We continue this procedure until we reach
level 1, moving back a level if ever the inequality in \eqref{eq:searchLevelK} for the current level is no longer satisfied. When we finish with level $1$, we will have found an integer point $\hat{x}$. Since $\hat{x}$ must satisfy \eqref{eq:searchIneq1}, we can use it to define a new, smaller radius which will effectively shrink the search space since there will be fewer valid choices for values of $x$.  To do this, update $\beta$ as follows, $$\beta = \left \| \bar{y} - R\hat{x} \right \|_2$$ and try to find a better integer point which lies in the new, smaller ellipsoid. We do this by continuing the search process from the current level after the radius has been tightened. However, at the current level $1$, there will be no integer which satisfies the inequality \eqref{eq:searchLevelK}, therefore we will have to start moving back up the search tree. Finally in the search process, when we can no longer find any $x_n$ to satisfy \eqref{eq:searchLevelK} with $k=n$, the search is complete and the last integer point $\hat{x}$ found is the solution.

\begin{figure}
\centering
\includegraphics[scale=0.5]{searchtree.png}
\caption{An example of the search process with solution $x = [-1,3,1]^T$.}
\label{fig:treeSearch}
\end{figure}

The above search process is actually a depth-first tree search in a tree of height $n$, see Figure \ref{fig:treeSearch},
where the number in a node  denotes the step number at which the node is encountered. Each edge going from level $i$ to $i-1$ represents fixing $x_{i-1}$ to some value, and each edge has a weight which is given by $r_{ii}^2(z_i -
c_i)^2$. We call the sum of the edge weights from the root to any node in the tree the ``cost'' to visit that node. This cost should not be confused with computational cost, when we discuss computational cost it will be mentioned explicitly. Notice that if we take the sum of these weights from $i$ to $n$ we simply get the partial
residual squared for fixing $x_{i:n}$, which is defined as $$\left \| \bar{y}_{i:n} - R_{i:n,i:n}x_{i:n} \right \|_2^2.$$ 

In this tree we define all leaves to have depth $n$. For example, we would not consider node $6$ in Figure \ref{fig:treeSearch} a leaf, it does not have children in the figure simply because none of its possible children satisfied the constraint from the search process \eqref{eq:searchIneq1}. This implies that each leaf in the tree
represents an integer vector $x \in \mathbb{Z}^n$ and its cost is the squared residual $$\left \|
 \bar{y} - Rx \right \|_2^2.$$ This means that the OILS problem is equivalent to
finding the lowest cost leaf in the tree. The SE enumeration takes advantage
of the fact that we can easily calculate the lowest cost child of any given
node to make the depth first search more efficient. The ILS problem is different from a classical tree search problem where in order to find the lowest weight edge leaving a given node you have to check each edge individually. In fact we can easily visit
the children of a node in order of increasing cost, that is what we are doing when we initially choose $x_k = \lfloor c_k \rceil$, and next choose it to be the second and third nearest integer to $c_k$.

A modification of the SE enumeration can be used to solve the BILS problem. To ensure that we remain within the box constraint, instead of choosing $x_k = \lfloor c_k \rceil$ at step $k$, we choose $x_{k} = \lfloor c_{k} \rceil_{{\cal B}_{k}}$, where ${\cal B}_{k}$ comes from \eqref{eq:boxCon}. Also recall that the notation $\lfloor c_{k} \rceil_{{\cal B}_{k}}$ denotes rounding $c_k$ to the nearest integer in the set ${\cal B}_{k}$. Suppose $x_{k-1:n}$ are fixed, then we must also ensure that as we visit the node corresponding to the second nearest integer (and all subsequent integers) for $x_k$ that we remain within the box constraint. This is trivial to accomplish, we simply stop incrementing $x_k$ if we hit the upper bound, and stop decrementing it if we hit the lower bound. If all values for $x_k$ that are within the box constraint have been used up but we are still within the area defined by the ellipsoid, we move back to level $k-1$. Following this process will yield the optimal BILS solution. For more information on the BILS implementation of the SE algorithm, see \cite{ChaH05}.

The choice of the initial value for $\beta$ can have a large impact on the number of nodes visited during the search. If $\beta$ is initially chosen to be infinity, the first point found by the SE search is known as the Babai integer point \cite{Bab86}. In some literature the Babai point refers to the corresponding lattice point $Rx$. One simple method for choosing an initial $\beta$ which may provide better results than infinity is as follows. Let $x_{\scriptscriptstyle{\mathrm{RLS}}}$ be the ordinary real LS solution. Then we may use the following initial $\beta$ for BILS problems,  $$\left \| R \lceil x_{\scriptscriptstyle{\mathrm{RLS}}} \rfloor_{{\cal B}} - \bar{y} \right \|_2,$$ for OILS problems we may use  $$\left \| R \lceil x_{\scriptscriptstyle{\mathrm{RLS}}} \rfloor - \bar{y} \right \|_2.$$ Obviously the optimal solution must have a residual less than or equal to these initial values for $\beta$. While the search process is never slower when using this initial $\beta$, often it will just end up being equivalent to using infinity since the Babai integer point is usually a better estimate than the rounded real least squares solution. In Chapter \ref{chap:Reduction} another idea for choosing an initial $\beta$ will be introduced.

\chapter{Reduction Strategies} \label{chap:Reduction}

This chapter focuses on the reduction stage for solving ILS problems. The goal is to modify the matrix $H$ and vector $y$ in such a way to preserve the original solution but make the problem easier. First the types of operations that we may perform on $H$ and $y$ will be described. The operations that are allowed are different for OILS and BILS problems, the difference and reason for it will be explained. Next some BILS reduction strategies will be described, we will show that two of the most effective strategies in the literature are actually equivalent and present a new algorithm that is based on combining ideas from these two which is more efficient than the originals and numerically reliable. Finally, it was found that some of the ideas from the BILS reduction algorithms can be applied to the OILS problems as well. This idea will be explored and some numerical results presented.

Consider the OILS problem defined in equation \eqref{eq:ils0}. The goal of the reduction is to modify the matrix $H$ in such a way that we still obtain the same solution $x$, but in fewer steps. With this goal, it is essential to know what types of operations we are allowed to perform on the matrix $H$ so that the solution $x$ is not modified. The first type of operation to consider is the orthogonal transformation. Suppose we apply some orthogonal matrix $Q$ to $H$ and $y$ from the left, then it is easy to see that the value of the objective function in equation \eqref{eq:ils0} remains unchanged. 

The second type of transformation that we may apply to the matrix $H$ from the right is any unimodular matrix. Unimodular matrices are square, integer matrices with determinant $1$ or $-1$. We can show that the inverse of a unimodular matrix is an integer matrix through basic linear algebra; this property is crucial for our application. Consider applying such a unimodular matrix $Z$ to $H$ from the right, equation \eqref{eq:ils0} becomes: 
$$
\| y - Hx  \|_2 = \| y - HZZ^{-1}x  \|_2 = \| y - \bar{H}z \|_2,
$$
where $\bar{H} = HZ$ and $z = Z^{-1}x$. When we solve this new OILS problem, we obtain some solution $z$. If $Z$ is a known unimodular matrix, we can solve for $x$ by computing $x = Zz$, this gives us the OILS solution to the original problem. Note that the unimodular matrix $Z$ guarantees the transformation from integer vectors to integer vectors.

Finally, it is worth mentioning that permutation matrices are unimodular since they will be used extensively throughout this chapter. It is obvious that the effect on the solution vector $x$ from applying a permutation matrix to $H$ from the right is to reorder the elements of the the solution $x$, this is important when considering the BILS problem.

When reducing the BILS problem we must be more careful. Consider the constraints on the solution $x$, defined in equation \eqref{eq:boxCon}. Applying orthogonal matrices to $H$ from the left has no effect on $x$, so the constraints are also unaffected. Applying permutation matrices to $H$ from the right will reorder the elements of the solution $x$, so we must reorder the elements in the constraint vectors as well, which is a trivial operation. However if we apply a general unimodular matrix $Z$ to $H$ from the right, the simple box constraints become complicated; we can no longer easily enforce the constraints on the solution during the search process. To see why this is the case, recall the SE search process. Suppose the search is currently at level $n$ and we would like to know to which value we should fix $z_{n}$ (note that $z$ denotes $Z^{-1}x$). We know that $$l_{n} \le x_{n} \le u_{n},$$ but the bound for $z_{n}$ is more complicated, $$l_{n} \le Z_{n,:}z \le u_{n}.$$ We would need to know the entire vector $z$ in order to compute $Z_{n,:}z$. At this point in the search, we have not yet fixed any of the values in $z$. The only way to complete the search process is to ignore the bounds on $z$, find a potential solution, compute $x = Zz$, and then check if each element of $x$ is within the bounds. This is extremely inefficient since many potential solutions could be eliminated very early on in the search if we were able to enforce the constraints during the search process. It is for this reason that when reducing BILS problems, we only consider orthogonal transformations and column permutations.

\vsp \section{BILS Reduction Algorithms} \label{sec:BILSReduction}

Due to the difficulty of applying general unimodular matrices to reduce the BILS problem, the algorithms in this section will focus on finding some permutation of the columns of the matrix $H$ in order to optimize the search process.

\vsp \subsection{Previous Reductions}

The ``Vertical-Bell Laboratories Layered Space-Time'' or V-BLAST algorithm \cite{FosGVW99} mentioned in section \ref{sec:prevWork} is a commonly used strategy to calculate the column permutations for $H$. Suppose we are working with $R$, which comes from the QR decomposition of $H$. Recall that the product of the diagonal elements of an upper triangular matrix is equal to the its determinant, and this value is invariant under permutation. This means that any time we swap two columns in the matrix $H$ and re-calculate the QR decomposition, one diagonal element of $R$ will always increase and the other will decrease.

The goal of the V-BLAST algorithm is to proceed from $k=n, \dots,1$ and find the column $h_p$ from the set of columns $h_1, \dots, h_k$ such that when $h_p$ and $h_k$ are swapped, the magnitude of the diagonal element $|r_{kk}|$ is maximal among the $k$ choices. There is an efficient algorithm to compute such a column ordering which is described in \cite{ChaP07}.

In \cite{ChaH05}, the ``Sorted QR Decomposition'' or SQRD column reordering strategy originally presented in \cite{WubBRKK01} for the same purpose as V-BLAST, was proposed for the purpose of BILS reduction. In the SQRD algorithm, we perform the QR decomposition from columns $k = 1,\dots, n$, at each step choosing as the $k^{th}$ column the column from the set $k,\dots,n$ which gives the smallest magnitude for the diagonal element $|r_{kk}|$. This should yield large $|r_{kk}|$ toward the end of the matrix since the product of the diagonal elements is a constant. Note that both SQRD and V-BLAST only use the information in the matrix $H$.

In \cite{SuW05}, Su and Wassell considered the geometry of the BILS
problem for the case that $H$ is nonsingular and proposed a new column reordering algorithm (to be called
the SW algorithm from here on for convenience) which uses all information of the BILS problem.
Unfortunately, the geometric interpretation of this algorithm is difficult to understand.
Probably due to page limit, the description of the algorithm is very concise, 
making efficient implementation difficult for ordinary users.

This thesis will give some new insight of the SW algorithm from an algebraic point of view.
Some modifications will be made so that the algorithm becomes more efficient
and easier to understand and furthermore it can handle a general full column rank $H$.
It is worth mentioning that the SW algorithm is not numerically stable. The numerical stability is not
necessarily crucial since a wrong answer just results in a different set of
permutations for the columns of H where any set of permutations is allowable.
Until this point, other column reordering algorithms only considered the matrix $H$.

Independently  Chang and Han in \cite{ChaH05} proposed
another column reordering algorithm (which will be referred to as  CH).
Their algorithm also uses all information of the BILS problem and the derivation
is based on an algebraic point of view. It is  easy to see from the equations in
the search process exactly what the CH column reordering is doing and why we
should expect a reduced complexity in the search process. The detailed
description of the CH column reordering is given in \cite{ChaH05} and it is easy
for others to implement the algorithm.
But our numerical tests and theoretical analysis indicated that CH has a higher complexity than SW, when SW
is implemented efficiently.
Our numerical tests also showed that CH and SW {\em almost} always   
produced the same permutation matrix $P$, the only exception is when then matrix $H$ is very ill conditioned.

In this section it will be shown that the CH algorithm and the  (modified)  SW algorithm give the same
column reordering in theory; in practice the reordering can be different due to the numerical instability of the SW algorithm. This is interesting because both algorithms were derived through different motivations
and we now have both a geometric justification and an algebraic justification 
for why the column reordering strategy should reduce the complexity of the search.
Furthermore, using the knowledge that certain steps in each algorithm are equivalent,
we can combine the best parts from each into a new algorithm. The new algorithm
has a lower flop count than either of the originals.
This is important to the successive interference cancellation decoder, 
which computes a suboptimal solution to the BILS problem. When computing suboptimal solutions, the overall runtime may not be so dominated by the search, and the complexity of the reduction can become a bottleneck.
The new algorithm can be interpreted in the same way as CH,
so it is easy to understand. 

In the following subsections, the CH and SW algorithms will be described in detail. This is necessary to understanding the proof of their equivalence and to see the motivation for the new algorithm. Also, a new algebraic interpretation and some improvements will be given for the SW algorithm. Finally the proof of equivalence of CH and SW will be given, and the new algorithm will be presented.

\vsp \subsection{CH Algorithm} \label{subsec:CH}
The CH algorithm first computes the QR decomposition of $H$,
then  tries to reorder the columns of $R$.
The motivation for this algorithm comes from observing equation \eqref{eq:searchLevelK}.
If the inequality is false we know that the current choice for the value of
$x_k$ given $x_{k+1:n}$ are fixed is incorrect and we prune the search tree (we do not need to explore the subtree for the current node). We
would like to choose the column permutations so that it is likely that the
inequality will be false at higher levels in the search tree, this way we waste less time exploring solutions that are not optimal. The CH column reordering strategy
does this by trying to maximize the left hand side of  \eqref{eq:searchLevelK} with large values of $\left | r_{kk}
\right |$ and minimize the right hand side by making $\left | r_{kk}(x_k-c_k) \right |$
large for values of $k = n,n-1, \dots, 1$.

Here we describe  step 1 of the CH algorithm, which determines the last column of the final $R$ 
(or equivalently the last column of the final $H$).
Subsequent steps are the same but are applied to a subproblem that is one dimension smaller. 
In step 1, for $i = 1,\dots,n$ we interchange
columns $i$ and $n$ of  $R$ (thus entries of $i$ and $n$ in $x$ are also swapped), then return $R$ to upper-triangular
by a series of Givens rotations applied to $R$ from the left, which  are also applied to $\bar{y}$. The following example demonstrates the process of returning the matrix $R$ to upper triangular after column $5$ is swapped with column $2$ in a $5 \times 5$ matrix.

\begin{equation} \label{eq:swappedCols}
\begin{bmatrix}
\times & \times & \times & \times & \times\\ 
  & \times & \times & \times & \times\\ 
  & \times & \times & \times &  \\ 
  & \times &  &  \times &  \\ 
  & \times &  &    & 
\end{bmatrix}
\end{equation}

Equation \ref{eq:swappedCols} shows the matrix directly after the column swap. We want to restore it to upper triangular. To do so, we start by using $n-i$ Givens rotations to zero the subdiagonal elements in column $i$, in this case $i = 2$. Each Givens rotation is an orthogonal matrix which adds multiples of two rows to eachother. For our purposes, we would like to add only multiples of adjacent rows. A Givens rotation that uses row $k$ to zero element $j$ in row $k+1$ is defined in equation \eqref{eq:givensRotation}, denote such a Givens rotation as $G_{k,k+1}$:

\begin{equation} \label{eq:givensRotation}
\begin{bmatrix}
I_{k-1} &  &  & \\ 
 & c & s & \\ 
 & -s & c & \\ 
 &  &  & I_{n-k-1}
\end{bmatrix} ,\quad c^2+s^2=1
\end{equation}
where $c=\frac{r_{k,j}}{\sqrt{r_{k,j}^2+r_{k+1,j}^2}}$ and $s=\frac{r_{k+1,j}}{\sqrt{r_{k,j}^2+r_{k+1,j}^2}}.$

We can use rotations $G_{n-1,n}, G_{n-2,n-1}, \dots, G_{i,i+1}$ to zero the subdiagonal elements in the $i^{th}$ column, however this creates subdiagonal elements in columns $i+1, \dots, n$. Equation \eqref{eq:subDiagEntries} shows the matrix after applying this first round of Givens rotatations.

\begin{equation} \label{eq:subDiagEntries}
\begin{bmatrix}
\times & \times & \times & \times & \times \\ 
  & \times & \times & \times & \times \\ 
  &   & \times & \times & \times \\ 
  &   & \times & \times &   \\ 
  &   &  &  \times & 
\end{bmatrix}
\end{equation}

We can now use a second round of Givens rotation to eliminate the new subdiagonal entries and restore the matrix to upper triangular. We use the rotations in the following order, $$G_{i+1,i+2}, G_{i+2,i+3}, \dots, G_{n-1,n}$$ where each rotation should zero the subdiagonal element in the corresponding column.

To avoid confusion, we denote the new $R$ after restoring the upper triangular structure by $\hat{R}$ and the new $\bar{y}$ by $\hat{y}$, where $\hat{y}$ is $\bar{y}$ with the Givens rotations applied to it.
We then compute  $c_n=\hat{y}_n/\hat{r}_{n,n}$ and 
\be
x_i^c=\arg\min_{x_i\in {\cal B}_i}\left |\hat{r}_{nn}(x_i- c_n) \right | = \left \lfloor c_n \right \rceil_{{\cal B}_i},
\label{eq:xic}
\ee
where the superscript $c$ denotes the CH algorithm. 
Let $\bar{x}_i^c$ be the second closest integer in ${{\cal B}_i}$ to $c_n$,
i.e.,  
$$\bar{x}_i^c= \left \lfloor c_n \right \rceil_{{\cal B}_i\backslash x_i^c}.$$
%So $\bar{x}_i^c = x_i^c \pm 1$.
Define
\be
\dist_i^c = |\hat{r}_{nn}( \bar{x}_i^c -c_n) |, 
\label{eq:dic}
\ee
which represents the partial residual given when $x_i$ is taken to be $\bar{x}_i^c$.
Let $j = {\arg\max}_i \dist_i^c$.
Then  column $j$ of the original $R$ is chosen to be the $n^{th}$ column of the final $R$.
%The algorithm then applies the Givens rotations that were used to restore $R$ to
%upper triangular to $\bar{y}$ and 
With the corresponding updated upper triangular $R$ and $\bar{y}$
(here for convenience we have removed hats),
the algorithm then updates  $\bar{y}_{1:n-1}$ again
by setting $$\bar{y}_{1:n-1}: = \bar{y}_{1:n-1} - r_{1:n-1,n}x_j$$ where $x_j=x_j^c$. 
%Note that we use $z$ to denote the permuted $x$.
Choosing  $x_j$ to be  $x_j^c$ here is exactly the same as what the search process does.
We then continue to work on the subproblem 
\be
\min_{\tilde{x}\in \mathbb{Z}^{n-1}} \left \| \bar{y}_{1:n-1}-R_{1:n-1,1:n-1}\tilde{x} \right \|_2,
\label{eq:subc}
\ee
where $$\tilde{x}=[x_1,\ldots, x_{j-1}, x_n, x_{j+1}, \ldots x_{n-1}]^T$$ satisfies the corresponding box constraint.
The pseudocode of the CH algorithm is given in Algorithm \ref{alg:CH}.

To determine the last column, CH finds the permutation to 
maximize $\left |r_{nn}(\bar{x}_i^c-c_n) \right |$. Using $\bar{x}_i^c$ instead of $x_i^c$
ensures that $\left | \bar{x}_i^c-c_n \right |$ is never less than $0.5$ but
also not very large. This means that usually if $\left | r_{nn}(\bar{x}_i^c-c_n)
\right |$ is large, $\left | r_{nn} \right |$ is large as well and the
requirement to have large $|r_{nn}|$ %in order from $m,\dots,1$  
is met.
Using $x_i^c$ would not be a good choice because $\left | x_i^c - c_n \right |$ might be 
very small or even $0$, then column $i$ would not be chosen to be column $n$
even if the corresponding $|r_{nn}|$ is large and on the contrary a column with small $|r_{nn}|$
but large $|x_i^c-c_n|$ may be chosen.

Now we will consider the complexity of CH. The 
significant computational cost comes from line \ref{l:chg} in Algorithm \ref{alg:CH},
which requires $6(k-i)^2$ flops.
If we sum this cost over all loop iterations and add the cost of the QR decomposition by Householder transformations, 
we  get a total complexity of $0.5n^4+2mn^2$ flops.

\begin{algorithm}
\caption{CH Algorithm - Returns $p$, the column permutation vector}
\label{alg:CH}
\begin{algorithmic}[1]
\STATE $p := 1:n$
\STATE $p' := 1:n$
\STATE Compute the QR decomposition of $H$: $\bsmx Q_1^T \\ Q_2^T \esmx H= \bsmx R\\ 0 \esmx$
             and compute  $\bar{y} : = Q_1^Ty$
  \FOR{$k=n$ to $2$}
  	\STATE $maxDist := -1$
    \FOR{$i=1$ to $k$}
    	\STATE $\hat{y} := \bar{y}_{1:k}$
    	\STATE $\hat{R} := R_{1:k,1:k}$
        \STATE  \label{l:chg} swap columns $i$ and $k$ of $\hat{R}$, return it  to upper
triangular with Givens rotations, also apply the Givens rotations to $\hat{y}$ %\hfill ($6(k-j)^2$ flops)  
        \STATE $x_i^c := \left \lfloor\hat{y}_k/\hat{r}_{k,k}\right
\rceil_{{\cal B}_i}$
        \STATE $\bar{x}_i^c := \left \lfloor\hat{y}_k/\hat{r}_{k,k}\right
\rceil_{{\cal B}_i\backslash x_i^c}$
        \STATE $dist_i^c := \left | \hat{r}_{k,k}\bar{x}_i^c - \hat{y}_k
\right | $
        \IF{$dist_i^c > maxDist$}
        	\STATE $maxDist := dist_i^c$
        	\STATE $j := i$
        	\STATE $R' := \hat{R}$
        	\STATE $y' := \hat{y}$
        \ENDIF
    \ENDFOR
    \STATE $p_k := p'_j$
    \STATE Interchange the intervals ${{\cal B}_k}$ and ${{\cal B}_{j}}$
    \STATE Intechange entries $k$ and $j$ in $p'$
    \STATE $R_{1:k,1:k} := R'$
    \STATE $\bar{y}_{1:k} := y' - R'_{1:k,k}x_j^c$
  \ENDFOR
  \STATE $p_1 := p'_1$
\end{algorithmic}
\end{algorithm}

\vsp \subsection{SW Original Algorithm} \label{subsec:SW}
The motivation for the SW algorithm comes from examining the geometry of the search process.

\ifx\du\undefined
  \newlength{\du}
\fi
\begin{figure}
\centering
\includegraphics[scale=0.7]{sumotivation.png}
\caption{Geometry of the search  with two different column ordering.}
\label{SEGeometry}
\end{figure}

Figure \ref{SEGeometry} shows a 2-D BILS problem;
\ref{SEGeometry}(a) represents the original column ordering and 
\ref{SEGeometry}(b) is after the columns have been swapped.

In the SW algorithm $H=[h_1,\ldots, h_n]$ is assumed to be square and non-singular.
Let 
$$
G =[g_1,\ldots, g_n]= H^{-T}.
$$ 
For any  integer $\alpha$,   \cite{SuW05} defines the
affine sets, $$F_i(\alpha) = \{w \ | \ g_i^T(w-h_i\alpha) = 0\}.$$
%which can also be thought of as $h_i$ shifted by the integer $\alpha$.
The lattice points generated by $H$ occur at the intersections of these affine sets. 
Let the orthogonal projection of a vector $s$ onto a vector $t$ be denoted as
$\mbox{proj}_t(s)$, then %\cite{SuW05} claims that 
the orthogonal projection of some vector $s$ onto $F_i(\alpha)$ is 
$$\mbox{proj}_{F_i(\alpha)}(s) = s -\mbox{proj}_{g_i}(s-h_i\alpha).$$ 
Therefore the orthogonal distance between $s$ and
$F_i(\alpha)$ is $$\dist(s,F_i(\alpha)) =  \| s - \mbox{proj}_{F_i(\alpha)}(s) \|_2.$$
%%%%%%%%%%%%
%\footnote{XW: Why don't they simply say the distance is just $\|\mbox{proj}_{g_i}(s-h_i\alpha)\|$?
%It's quite weird. SB: I'm not sure, it seems like that projection is equivalent
%and much simpler and we could avoid using the affine sets $F_i(s)$}
%%%%%%%%%%%%%%%%
In \cite{SuW05}, the points labeled $\mbox{proj}_{F_2(1)}(y)$ and
$\mbox{proj}_{F_2(-1)}(y)$ in Figure \ref{SEGeometry}
are called residual targets  and ``represent the
components [of $y$] that remain after an orthogonal part has been projected
away.''  

Note that $F_2(\alpha)$ in Figure \ref{SEGeometry} is a sublattice of dimension $1$.  
Algebraically it is the lattice generated by $H$ with column
$2$ removed. It can also be thought of as a subtree of the search tree where
$x_2 = \alpha$ has been fixed. 
In the first step of the search process for a general case, $x_n$ is chosen to be
$$x_n=\arg\min_{\alpha\in {\cal B}_n}\dist(y,F_n(\alpha));$$ thus $F_n(x_n)$ is the nearest affine set to $y$. 
Actually the value of $x_n$ is identical to $\lfloor c_n \rceil_{{\cal B}_n}$ given in Chapter ~\ref{chap:SESearch},
which will be proven later.
Then  $y$ is updated  as $$y := \mbox{proj}_{F_n(x_n)}(y) - h_nx_n.$$ If we look
at Figure \ref{SEGeometry}, we see that the projection $\mbox{proj}_{F_n(x_n)}(y)$ moves $y$ onto $F_n(x_n)$, while the subtraction of $h_nx_n$ algebraically fixes the value of $x_n$. This is necessary because in subsequent steps we will not consider the column $h_n$.

We now apply the same process to the new $n-1$ dimensional search space
$F_n(x_n)$. If at some level $i$, $$\min_{\alpha\in {\cal B}_i}\dist(y,F_i(\alpha))$$ exceeds the current
search radius, we must move back to level $i+1$. % and choose $x_i = x_i \pm 1$.
When the search process reaches level $1$ and fixes $x_1$, it updates the radius to  
$\dist(y,F_1(x_1))$ and moves back up to level $2$.

Note that this search process is mathematically equivalent to the one described in Chapter
\ref{chap:SESearch}; the difference is that it  does projections
because  the generator matrix is not assumed to be upper-triangular. 
Computationally the former is more expensive than the latter since the cost for QR decomposition is only incurred once and we avoid the expensive projections.

To see the motivation of the SW algorithm for choosing a particular column ordering,
consider Figure \ref{SEGeometry}. Suppose the search algorithm has knowledge of
the residual for the optimal solution (the radius of the circle in the diagram).
With the column ordering chosen in (a), there are two possible choices for $x_2$,
leading to the two dashed lines $F_2(-1)$ and $F_2(1)$ which cross the circle. This means
that we will need to compute $x_1$ for both of these choices
before we can determine which one leads to the optimum solution. In (b), there
is only one possible choice for $x_1$,  leading to the only dashed line $F_1(-1)$
which crosses the circle, meaning we only need to find $x_2$ to find the optimum solution.
Since the projection resulting from the correct choice of $x_2$ will always be
within the sphere, it makes sense to choose the ordering which maximizes the
distance to the second best choice for $x_2$ in hopes that the second nearest
choice will result in a value for $$\min_{\alpha\in {\cal B}_2}\dist(y,F_2(\alpha))$$ outside the sphere and the
dimensionality can be reduced by one. 
For more details on the geometry, see  \cite{SuW05}.

The following will give an overview of the SW algorithm as given in
\cite{SuW05} but described in a framework similar to what was used to describe
CH. In the first step to determine the last column, for each $i = 1, \dots, n$,
we compute 
\be
x_i^s \!=\! \arg\min_{\alpha\in {\cal B}_i}\dist(y,F_i(\alpha)) 
\!=\!  \arg\min _{\alpha\in {\cal B}_i}| y^Tg_i - \alpha | \!=\! \left \lfloor y^Tg_i \right \rceil_{{\cal B}_i},
\label{eq:xis}
\ee
where the superscript $s$  stands for the SW algorithm.
Let $\bar{x}_i^s$ be the second closest integer in ${{\cal B}_i }$ to $y^Tg_i$, i.e.,
$$\bar{x}_i^s = \left \lfloor y^Tg_i \right \rceil_{{\cal B}_i\backslash x_i^s}.$$
Let $j = \arg\max_i  \dist (y,F_i(\bar{x}_i^s))$. 
Then SW chooses column $j$ as the last column of the final reordered $H$,
updates $y$ by setting 
$$y:=\mbox{proj}_{F_j(x_j^s)}(y) -h_jx_j^s$$  and 
updates $G$ by setting $$g_i: = \mbox{proj}_{F_j(0)}(g_i)$$ for all $i\neq j$. 
After $G$ and $y$ have been updated, the algorithm continues to find column $n-1$ in the
same way etc. The pseudo-code of the SW algorithm is given in Algorithm \ref{alg:SWOrig}.
\begin{algorithm}
\caption{SW Algorithm - Returns $p$, the column permutation vector}
\label{alg:SWOrig}
\begin{algorithmic}[1]
\STATE $p := 1:n$
\STATE $p' := \{1, 2, \ldots, n\}$
\STATE \label{l:swG} $G := H^{-T}$ \hfill %($2n^3$ flops)
\FOR{$k=n$ to $2$}
	\STATE $maxDist := -1$
	\FOR{$i \in p'$}
		\STATE $x_i^s := \left \lfloor  y^Tg_i \right \rceil_{{\cal B }_i}$ \hfill% ($2n$ flops)
		     \label{l:swx}
		\STATE $\bar{x}_i^s := \left \lfloor y^Tg_i \right \rceil_{{{\cal B }_i}{\backslash x_i^s}}$
		    \label{l:swbx}
		\STATE  \label{l:swd} $\dist_i^s := \dist(y,F_i(\bar{x}_i^s))$ %\hfill ($3n^2$ flops)
		\IF{$dist_i^s > maxDist$}
			\STATE $maxDist := dist_i^s$
			\STATE $j := i$
		\ENDIF
	\ENDFOR
	\STATE $p_k := j$
	\STATE $p' := p' \backslash j$
	\STATE $y := \mbox{proj}_{F_j(x_j^s)}(y) -h_jx_j^s$  \label{l:swy}
	\FOR{$i \in p'$}
		\STATE \label{l:swg} $g_i := \mbox{proj}_{F_j(0)}(g_i) $ %\hfill ($2n^2$ flops) 
	\ENDFOR
\ENDFOR
\STATE $p_1 := p'$
\end{algorithmic}
\end{algorithm}

%Now we look at the complexity of Algorithm \ref{alg:SWOrig}.
Su and Wassell did not say how to implement the algorithm and did not give a complexity analysis.
The parts of the computational cost we must consider for implementation occur in
lines \ref{l:swd} and \ref{l:swg}.
Note that $$\dist(y,F_i(\bar{x}_i^s))=\|\mbox{proj}_{g_i}(y-h_i\bar{x}_i^s)\|_2$$
and $$\mbox{proj}_{F_j(0)}(g_i)= g_i -\mbox{proj}_{g_i} g_i,$$
where $$\mbox{proj}_{g_i}=g_ig_i^\dag = g_ig_i^T/\|g_i\|^2.$$ 
A naive implementation would first compute $\mbox{proj}_{g_i}$, requiring $n^2$ flops, then compute 
$$\|\mbox{proj}_{g_i}(y-h_i\bar{x}_i^s)\|_2$$ and $$g_i -\mbox{proj}_{g_i} g_i,$$ each requiring $2n^2$ flops.
Summing these costs over all loop iterations we get a total complexity of $2.5n^4$ flops.
In the next subsection we will simplify some steps in Algorithm \ref{alg:SWOrig}
and show how to  implement them efficiently.

\vsp \subsection{SW Algorithm Interpretation and Improvements}
\label{sec:improvedSW}
In this section we give a new algebraic interpretation of some steps in Algorithm 2,
simplify some key steps to improve the efficiency,
and  extend the algorithm to handle a more general case.
All line numbers refer to Algorithm \ref{alg:SWOrig}.

%At each step, we must compute  $x_i^s = \left \lfloor y^Tg_i \right \rceil_{{\cal B}_i }$ and $\dist_i^s$. 
First we show  how to efficiently compute $\dist_i^s$ in line \ref{l:swd}. 
Observing that $g_i^Th_i = 1$, we have
\begin{equation}
\label{eq:newDist}
\dist_i^s =   \|  g_ig_i^\dag (y-h_i\bar{x}_i^s)  \|_2 
=    | y^Tg_i -\bar{x}_i^s |/\| g_i  \|_2.
\end{equation} 
Note that $y^Tg_i$ and $\bar{x}_i^s$ have been computed in lines \ref{l:swx} and \ref{l:swbx}, respectively.
So the main computational cost of computing $\dist_i^s$ is the cost of computing $\|g_i\|_2$,
requiring only $2n$ flops. 
For $k=n$ in Algorithm 2,  $$y^Tg_i=y^TH^{-T}e_i=(H^{-1}y)^Te_i,$$ i.e.,  $y^Tg_i$
is the $i^{th}$ entry of the real solution $x$ for $Hx=y$. 
The interpretation can be generalized to  an arbitrary $k$. 

In line \ref{l:swg} of Algorithm 2,  
\begin{align}
g_i^{\small \mbox{new}} & \equiv \mbox{proj}_{F_j(0)}(g_i)   \nonumber \\
  & =(I- \mbox{proj}_{g_j})g_i=g_i- g_j(g_j^Tg_i/\|g_j\|_2^2). \label{eq:gup}
\end{align}
Using the last expression for computation needs only $4n$ flops
(note that $\|g_j\|_2$ has been computed before, see \eqref{eq:newDist}).
We can actually show that the above is performing updating of $G$, the Moore-Penrose generalized inverse of
$H$ after we remove its $j^{th}$ column. This will allow us to interpret the subsequent steps of SW as operating on a subproblem. For proof of this, see \cite{Cli64}.

In line \ref{l:swy} of Algorithm 2,
\begin{align}
y^{\small \mbox{new}} & \!\equiv\! \mbox{proj}_{F_j(x_j^s)}(y) - h_jx_j^s 
 \!=\! (y -  g_jg_j^\dag(y-h_jx_j^s)) - h_jx_j^s  \nonumber \\
   &=  (I-\mbox{proj}_{g_j})(y-h_jx_j^s). \label{eq:yup}  
\end{align}
This means that after $x_j$ is fixed to be $x_j^s$, $h_jx_j^s$ is combined with $y$ (the same
as CH does)  and then the vector is projected to the orthogonal complement of 
the space spanned by $g_j$. 
We can show that this guarantees that the updated $y$ is in the subspace spanned by
the columns of $H$ which have not yet been chosen.
This is consistent with the assumption that $H$ is nonsingular, which implies that 
the original $y$  is in the space spanned by  the columns of $H$.
However, it is not necessary to apply the orthogonal projector $I- \mbox{proj}_{g_j}$ to $y-h_jx_j^s$ in \eqref{eq:yup}.
The reason is as follows. 
In Algorithm 2, $y^{\small \mbox{new}}$ and $g_i^{\small \mbox{new}}$ will be used only for computing 
$(y^{\small \mbox{new}})^Tg_i^{\small \mbox{new}}$ (see line \ref{l:swx}).
But from \eqref{eq:gup} and \eqref{eq:yup}
\begin{align*}
(y^{\small \mbox{new}})^Tg_i^{\small \mbox{new}}
& =(y-h_jx_j^s)^T(I-\mbox{proj}_{g_j})(I-\mbox{proj}_{g_j})g_i \\
&=(y-h_jx_j^s)^Tg_i^{\small \mbox{new}}.
\end{align*}
Therefore, line \ref{l:swy} can be replaced by $y:=y-h_jx_j^s$.
This not only simplifies the computation but also is much easier to interpret---after $x_j$ is fixed to be $x_j^s$,  
$h_jx_j^s$ is combined into $y$ as in the CH algorithm.
Let $H_{:,1:n-1}$ denote $H$ after its $j^{th}$ column is removed. 
We then continue to work on the subproblem
\be
\min_{\check{x}\in \mathbb{Z}^{n-1}}\|y-H_{:,1:n-1}\check{x}\|_2, 
\label{eq:subs}
\ee
where $$\check{x}=[x_1,\ldots,x_{j-1},x_{j+1},\ldots,x_n]^T$$ 
satisfies the corresponding box constraint.
Here $H_{:,1:n-1}$ is not square. But there is no problem to handle
it, as explained in the next paragraph.

In \cite{SuW05},  $H$ is assumed to be square and non-singular. 
In our opinion, this condition may cause confusion,
since for each $k$ except $k=n$ in Algorithm 2, 
the remaining columns of $H$ which have not been chosen do not form a square matrix.
Also the condition restricts the application of the algorithm to a general full column rank matrix $H$,
unless we transform $H$ to a nonsingular matrix $R$ by the QR decomposition.
To extend the algorithm to a general full column rank matrix $H$, we need only 
replace line \ref{l:swG} by $G:=(H^{\dagger})^T$.
This extension has another benefit. 
We mentioned before that the updating of $G$ in line \ref{l:swg}
is actually the updating of the Moore-Pernrose generalized inverse 
of the matrix formed by the columns of $H$ which have not been chosen. 
So the extension makes all steps consistent.

To reliably compute $G$ for a general full column rank $H$,
we can compute the  QR decomposition $H=Q_1R$ by the Householder transformations
and then solve the triangular system  $RG^T=Q_1^T$ to obtain $G$.
This requires $(5m-4n/3)n^2$ flops. 
Another less reliable but more efficient way to do this is to compute $G=H(H^TH)^{-1}$. 
To do this efficiently we would compute the Cholesky decomposition  $H^TH = R^TR$ and solve 
$R^TRG^T = H^T$ for $G$ by using the triangular structure of $R$. 
The total computational cost for computing $G$ by this method can be shown to be $3mn^2+\frac{n^3}{3}$.
If $H$ is square and nonsingular, we would use the LU decomposition with partial pivoting to compute $H^{-1}$
and the cost is $2n^3$ flops.

For the rest of the algorithm if we use the simplification and efficient implementations
mentioned above, we can show that it needs $4mn^2$ flops. 

We see that the modified SW algorithm is much more efficient than both the CH algorithm
and the SW algorithm implemented in the naive way we mentioned in the previous subsection.

\vsp \subsection{Proof of Equivalence of SW and CH}
In this subsection we prove that  CH and  the modified  SW produce the same set of permutations
for a general full column rank $H$.
To prove this it will suffice to prove that $x_i^s = x_i^c$, $\bar{x}_i^s =\bar{x}_i^c$,
$\dist_i^s = \dist_i^c$ for $i=1, \ldots, n$ in the first step which determines the last column of the final reordered $H$
and that the subproblems produced for the second step of
each algorithm are equivalent. 

Proving $x_i^s = x_i^c$ is not difficult.
The only effect the interchange of columns $i$  and $n$ of $R$ in CH  
has on the real LS solution is that elements $i$ and $n$ of the solution are swapped.
Therefore $x_i^c$ is just the $i^{th}$ element of the real LS
solution rounded to the nearest integer in ${{\cal B}_i }$. 
Thus, with equations \eqref{eq:xic} and \eqref{eq:xis},
\be
x_i^c=   \lfloor (H^{\dagger}y)_i  \rceil_{{\cal B}_i }
=  \lfloor e_i^T H^{\dagger}y   \rceil_{{\cal B}_i }
=  \lfloor g_i^T  y \rceil_{{\cal B}_i } =x_i^s.
\label{eq:xics}
\ee
Therefore we also have $\bar{x}_i^c=\bar{x}_i^s$.

In CH, after applying a permutation $P$ to swap columns $i$ and $n$ of $R$,  
we apply $V^T$, a product of the Givens rotations, to bring $R$ back to a new upper triangular
matrix, denoted by $\hat{R}$, and also apply $V^T$ to $\bar{y}$, 
leading to  $$\hat{y} = V^T\bar{y}.$$
Thus  $$\hat{R}=V^T RP$$ and $$\hat{y} = V^T\bar{y}=V^TQ_1^Ty.$$
Then $$H=Q_1R= Q_1V\hat{R}P^T,$$ $$H^\dag= P\hat{R}^{-1}V^TQ_1^T,$$ 
$$g_i=(H^\dag)^Te_i=Q_1V\hat{R}^{-T}P^Te_i=Q_1V\hat{R}^{-T}e_n,$$
and $$\|g_i\|_2=\|\hat{R}^{-T}e_n\|_2=1/|\hat{r}_{nn}|.$$
Therefore, with \eqref{eq:newDist} and \eqref{eq:dic}
\begin{align}
\dist_i^s
&=\frac{ | y^Tg_i - \bar{x}_i^s   |}{  \| g_i   \|_2} 
=|\hat{r}_{nn}||y^TQ_1V\hat{R}^{-T}e_n- \bar{x}_i^s  |  \label{eq:disc} \\
& = |\hat{r}_{nn}|| \hat{y}_n/\hat{r}_{nn} - \bar{x}_i^s | 
 = |\hat{r}_{nn}(c_n-\bar{x}_i^c)| =\dist_i^c.  \nonumber
\end{align}

Now we consider  the subproblem \eqref{eq:subc} in CH and the subproblem \eqref{eq:subs} in SW.
We can easily show that $R_{1:n-1,1:n-1}$ in  \eqref{eq:subc} is the $R$-factor of the QR decomposition
of $H_{:,1:n-1}P$, where $H_{:,1:n-1}$ is the matrix given in \eqref{eq:subs}
and $P$ is a permutation matrix such that $\check{x}=P\tilde{x}$,
and that $\bar{y}_{1:n-1}$ in  \eqref{eq:subc} is the multiplication of the transpose of 
the $Q_1$-factor of the QR decomposition of $H_{:,1:n-1}P$ and $y$ in \eqref{eq:subs}.
Thus the two subproblems are equivalent. To see proof of this, simply observe the following:
$$
HP = QRP = QV^TVRP
$$
where $V$ is the product of Givens rotations which restores R to an upper triangular matrix after the permutation $P$. Therefore, $VRP$ is upper triangular and obviously equivalent to the matrix used in the second step of CH by its definition; it is just a permutation applied to $R$ which is then restored to upper triangular through the product of Givens rotations $V$. The equivalence of $y$ in the second step should be simple to deduce.

\vsp \subsection{New Algorithm} \label{subsec:newReduction}
Now that we know the two algorithms are equivalent, we can take the best
parts from both and combine them to form a new algorithm. 
The main computational cost in CH is to interchange the columns of $R$ and return it to
upper-triangular form using Givens rotations. 
When we determine the $k^{th}$ column,  we must do this $k$ times. 
We can avoid all but one of these column interchanges by computing $x_i^c$, 
$\bar{x}_i^c$ and $\dist_i^c$ directly using the equations from SW. 

After the QR decomposition of $H$, we  solve the reduced BILS problem \eqref{eq:ils}.
We need only consider how to determine the last column of the final $R$.
Other columns can be determined similarly. 
Here we use the ideas from SW.
Let $G=R^{-T}$, which is lower triangular.
By \eqref{eq:xics}, we compute for $i=1,\ldots, n$
\begin{align*}
& x_i = \left \lfloor \bar{y}^TG_{:,i} \right \rceil_{{\cal B}_i} 
=\left \lfloor \bar{y}_{i:n}^T G_{i:n,i} \right \rceil_{{\cal B}_i}, \;
\bar{x}_i=\left \lfloor \bar{y}_{i:n}^T G_{i:n,i} \right \rceil_{{\cal B}_i\backslash x_i}, \\
& \dist_i = |\bar{y}_{i:n}^T G_{i:n,i}-\bar{x}_i|/\|G_{i:n,i}\|_2.
\end{align*}

Let $j=\arg\max_{i} \dist_i$. We  take a slightly different approach to
permuting the columns than was used in CH. Once $j$ is determined, we  
set $$\bar{y}_{1:n-1} := \bar{y}_{1:n-1} - r_{1:n-1,j}x_j.$$ Then we simply remove
the $j^{th}$ column from $R$, and restore it to upper triangular using
Givens rotations. We then apply the same Givens rotations to the new
$\bar{y}$. In addition, we must also update the inverse matrix $G$. This is
very easy, we can just remove the $j^{th}$ column of $G$ and apply the same
Givens rotations that were used to restore the upper triangular structure of
$R$. To see this is true notice that removing column $j$ of $R$ is
mathematically equivalent to rotating $j$ to the last column and shifting
columns $j, j+1, \ldots, n$ to the left one position, since we will only consider columns
$1, 2, \ldots, n-1$ in subsequent steps. Suppose $P$ is the permutation matrix which
permutes the columns as described, and $V^T$ is the product of Givens
rotations to restore $R$ to upper-triangular. Let $\hat{R} = V^TRP$ and set
$\hat{G} = \hat{R}^{-T}$. Then
$$
\hat{G} = (V^TRP)^{-T} = V^TR^{-T}P = V^TGP.
$$
This indicates that  the same $V$ and $P$, which are used to transform $R$ to  $\hat{R}$, 
also transform $G$ to $\hat{G}$.
Since $\hat{G}$ is lower triangular, it is easy to verify that
$\hat{G}_{1:n-1,1:n-1} = \hat{R}^{-T}_{1:n-1,1:n-1}$.
Both $\hat{R}_{1:n-1,1:n-1}$ and $\hat{G}_{1:n-1,1:n-1}$ will be used in the next step.

After this, as in the CH algorithm, we continue to work on the subproblem of size $n-1$. 
The advantages of using the ideas from CH are that we always have a lower triangular $G$
whose dimension is reduced by one at each step
and the updating of $G$ is numerically stable as we use orthogonal transformations. The lower triangular structure of $G$ and the upper triangular structure of $R$ allow us to save computation as the algorithm progresses.
The pseudocode of the new algorithm is given in Algorithm  \ref{alg:NEW}.

\begin{algorithm}
\caption{New algorithm}
\label{alg:NEW}
\begin{algorithmic}[1]
\STATE  Compute the QR decomposition of $H$ by Householder transformations: 
$\bsmx Q_1^T \\ Q_2^T \esmx H= \bsmx R\\ 0 \esmx$  \\
             and compute  $\bar{y} : = Q_1^Ty$ \hfill ($2(m-n/3)n^2$ flops)
\STATE $G := R^{-T}$ \hfill ($\frac{n^3}{3}$ flops)
\STATE $p := 1:n$
\STATE $p' := 1:n$
\FOR{$k=n$ to $2$}
	\STATE $maxDist := -1$
	\FOR{$i=1$ to $k$}
	         \STATE $\alpha=y_{i:k}^TG_{i:k,i}$
	         \STATE $x_i := \left \lfloor \alpha \right \rceil_{{\cal B}_i}$ \hfill ($2(k-i)$ flops)
	         \STATE $\bar{x}_i := \left \lfloor \alpha \right \rceil_{{{\cal B}_i}\backslash x_i}$
	         \STATE $\dist_i =|\alpha-\bar{x}_i|/ \| G_{i:k,i} \|_2$ \hfill ($2(k-i)$ flops)
			 \IF{$dist_i > maxDist$}
			 	\STATE $maxDist : = dist_i$
			 	\STATE $j:=i$
			 \ENDIF	
	\ENDFOR
	\STATE $p_k := p'_j$
	\STATE Interchange the intervals ${{\cal B}_k}$ and ${{\cal B}_j}$
	\STATE Interchange entries $k$ and $j$ in $p'$
	\STATE Set $\bar{y}:=\bar{y}_{1:k-1} - R_{1:k-1,j}x_j$	
	\STATE Remove column $j$ of $R$ and $G$, and return $R$ and $G$ to upper and lower triangular by Givens rotations, respectively, and then remove the last row of $R$ and $G$. The same Givens rotations are applied to $\bar{y}$. \\ \hfill ($6k(k-j)$ flops)
\ENDFOR
\STATE $p_1 = p'_1$
\end{algorithmic}
\end{algorithm}

Here we consider the complexity analysis of the new algorithm. 
If we sum the costs in Algorithm \ref{alg:NEW} over all loop iterations,
we get a total of $\frac{7n^3}{3} + 2mn^2$ flops in the worst case. 
The worst case is very unlikely to occur, it arises when $j=1$ each iteration of the outer loop. In the average case
however, $j$ is around $k/2$ and we get an average case complexity of $\frac{4n^3}{3} + 2mn^2$ flops.
In both cases, the complexity is less than the complexity of the modified SW algorithm.

\vsp \section{Reduction for OILS Problems} \label{sec:oilsreduction}
For the OILS problem, the most common reduction strategy is to apply the ``Lenstra-Lenstra-Lovász'' or LLL reduction \cite{LenLL82} to the matrix $H$. There are a few ways to describe the LLL reduction process and what it means for a matrix $H$ to be LLL reduced. In this thesis, we will look at the LLL algorithm as a matrix factorization, 
$$\begin{bmatrix}
Q_1^T\\ 
Q_2^T
\end{bmatrix}
HZ = 
\begin{bmatrix}
R\\ 
0
\end{bmatrix},$$
where $Q = [Q_1, Q_2] \in \mathbb{R}^{m \times m} $ is orthogonal, $R$ is upper triangular, and $Z \in \mathbb{Z}^{n \times n}$ is unimodular. After this special $QRZ$ decomposition (the $QRZ$ decomposition is not unique), the matrix $R$ is LLL reduced i.e.,

\begin{align} \label{eq:LLLConditions1}
&\left | r_{k-1,j} \right | \le \frac{1}{2} \left | r_{k-1,k-1} \right |, \quad
j = k:n, k=2:n \\
\label{eq:LLLConditions2}
&\delta r_{k-1,k-1}^2 \le r_{k-1,k}^2 + r_{k,k}^2
\end{align}

From \eqref{eq:LLLConditions1} which is known as the ``size reduction''
condition and \eqref{eq:LLLConditions2} which is called the ``Lovasz''
condition, we can easily obtain the following inequality:
\begin{equation} \label{eq:LLLdiagonal}
\left | r_{k-1,k-1} \right | \le \frac{2}{\sqrt{4\delta -1}}\left | r_{k,k} \right |.
\end{equation}

The parameter $\delta$ can take on values from the interval $(\frac{1}{4}, 1]$. Looking at
equation \eqref{eq:LLLdiagonal}, we can obtain some sense of why using a LLL reduced
matrix $R$ in the search process should yield a performance improvement. In many
applications, we use $\delta=1$, this is the maximum value for $\delta$ and
forces the most strict ordering of the diagonal, however as $\delta$ increases the LLL algorithm may take longer to converge. From here on when we discuss the LLL algorithm, we will assume $\delta = 1$.

We know from previous discussion
about the box constrained reductions that it is desirable to have large diagonal
elements, with the largest possible diagonal elements toward the end, $|r_{11}|
<, \dots,< |r_{nn}|$. Unfortunately this ordering can not always be obtained. The inequality implied by the LLL conditions, \eqref{eq:LLLdiagonal} does not impose this ordering strictly but it does ensure that it is not off by too much. In \cite{XieCB11} the authors demonstrate that some other properties of the LLL reduction also help with the search process. One of these relates to the distribution of the diagonal. They show that the minimum and maximum diagonal element should not be too far apart. This is consistent with the requirement for large diagonal elements since the product of the diagonal is constant. They also show that we should try and reduce the super diagonal entries of the matrix, which the LLL reduction accomplishes. 

These properties allows us to prune the search
tree at higher levels without wasting computational effort on suboptimal
solutions. Unfortunately these properties are not always possible to obtain.
However, the equation \eqref{eq:LLLdiagonal} does give us a guarantee about the
relative sizes of the diagonal elements and the LLL conditions give some bounds
on the size of the super diagonal.

\vsp \subsection{Computing the LLL Reduction}

This section will give details on how the LLL reduction described above can be
computed. It is interesting to note that the LLL reduction is not unique, other
methods to compute it exist and may yield different but still valid results.

In the following paragraphs the matrix operations that will be used by the LLL
reduction will be described. Once the operations which the reduction uses have
been presented, the LLL algorithm will be described and pseudocode given to
compute it.

\vsp \subsubsection{Integer Gauss Transformations.} \label{subsec:IGT}

One special type of unimodular matrix is an integer Gauss transformation (IGT), which can be defined as follows:
\begin{equation}
Z_{ij} = I-\mu e_ie_j^T, \quad \mu \in \mathbb{Z}.
\end{equation}

We would like to know how these transformations affect an upper triangular matrix $R$. Suppose we apply such an IGT to $R$ from the right, this gives:
\begin{equation}
\bar{R} = RZ_{ij} = R - \mu Re_ie_j^T.
\end{equation}
The overall effect of this transformation on the matrix R is that the $j^{th}$ column has some integer multiple of the $i^{th}$ column subtracted from it, therefore:
\begin{equation}
\bar{r}_{kj} = r_{kj} - \mu r_{ki}, \quad k=1 , \dots, i.
\end{equation}
If we take $\mu = \lfloor \frac{r_{ij}}{r_{ii}} \rceil$, it should be clear that
$|\bar{r}_{ij}| \le \frac{1}{2}|r_{ii}|$, so given a particular column $j$ in an
upper triangular matrix $R$, we should be able to use $j-1$ IGTs to satisfy the
size reduction condition \eqref{eq:LLLConditions1}.

\vsp \subsubsection{Permutations.} \label{subsec:Perm}
After doing IGTs, there is no guarantee that the Lovasz condition,
\eqref{eq:LLLConditions2} will be satisfied; often it is not. In this case, we
must permute the columns of $R$ in order for the condition to hold. If
$|r_{k-1,k-1}| > \sqrt{r^2_{k-1,k} + r^2_{k,k}}$, then we permute columns
$k$ and $k-1$. After performing the column permutation, $R$ will no longer be
upper triangular. To restore the upper triangular structure of $R$, we can apply
Givens rotations as we did in section \ref{sec:BILSReduction}. In this case
however, only one Givens rotation will be required to zero a single sub diagonal
element.

After the permutation, the Lovasz condition, \eqref{eq:LLLConditions2} will
hold. After performing this permutation, we also have the guarantee that value
of $|r_{kk}|$ will increase and $|r_{k-1,k-1}|$ will decrease, therefore the
resulting matrix will have something closer to an increasing diagonal.

By using the above permutations and IGTs together, we can devise an algorithm to
satisfy both of the LLL conditions \eqref{eq:LLLConditions1} and
\eqref{eq:LLLConditions2}. We will start by letting $H=[Q_1, Q_2] \bmx R \\ 0
\emx$ denote the $QR$ decomposition of the matrix $H$. We will work with the
columns of $R$ from right to left, starting with column $k=n$. The idea is to
move to the left so that at any step $k$, the columns $k+1:n$ satisfy the LLL
conditions. In the $k^{th}$ step, we start by using IGTs to make sure column $k$
satisfies the size reduction condition \eqref{eq:LLLConditions1}, $|r_{ik}| <
\frac{1}{2}|r_{ii}|, \quad i = k-1:-1:1$. If the Lovasz condition in
\eqref{eq:LLLConditions2} holds, we move to column $k-1$, otherwise column $k-1$
and $k$ are swapped with a column permutation and $R$ is brought back to upper
triangular as described in subsection \ref{subsec:Perm}. After applying a column
permutation, we must move back to column $k+1$ since it is possible that the
permutation we applied will cause the conditions on the previous column to no
longer be satisfied. When we reach column $1$, we know the matrix $R$ must be
LLL reduced. Algorithm \ref{alg:LLL} describes this process. As mentioned
previously in subsection \ref{subsec:Reductions}, in theory this algorithm terminates in a
finite number of iterations (although this has not been proven for fixed precision arithmetic). There are also some bounds on the complexity, for more
detail on the complexity and convergence, see \cite{WubSJM11}.

\begin{algorithm}
\caption{LLL Algorithm - Returns R the LLL reduced upper triangular matrix and Z a product of IGTs and permutations}
\label{alg:LLL}
\begin{algorithmic}[1]
\STATE Compute the QR decomposition of $H$: $\bsmx Q_1^T \\ Q_2^T \esmx H= \bsmx R\\ 0 \esmx$
\STATE $Z := I_n$
\STATE $k := n$
\WHILE{$k \ge 2$}
	\IF{$k > n$}
		\STATE $k := n$	
	\ENDIF
	\STATE Compute IGTs so that column $k$ satisfies $|r_{ik}| < \frac{1}{2}|r_{ii}|, \quad i = k-1:-1:1$, apply the IGTs 		to Z and R
	\IF{$r_{k,k}^2 + r_{k-1,k}^2 < r_{k-1,k-1}^2$}
		\STATE Swap columns $k$ and $k-1$ in R
		\STATE Swap columns $k$ and $k-1$ in Z
		\STATE $k := k+1$
	\ELSE
		\STATE $k := k-1$
	\ENDIF
\ENDWHILE
\end{algorithmic}
\end{algorithm}

\vsp \subsection{New Reduction Algorithm} \label{subsec:newOILS}
In section \ref{sec:BILSReduction}, the motivation for the permutation based reductions was given. While these algorithms make use of the box constraint, the original motivation for them does not rely on it. We may apply these permutation based strategies to unconstrained problems by simply setting the constraints to -infinity and +infinity.

Applying a permutation based reduction algorithm directly to the matrix $H$ however may yield results that are much worse than those given by LLL on average, this was confirmed through numerical experiments. There are a few reasons for why this is the case. First, since the permutation based reductions do not use IGTs, we can not reduce the super diagonal entries. In \cite{XieCB11} it was shown that reducing super diagonal entries is important for the search process. Also, we can not control the distribution of the diagonal entries without IGTs, and using only permutations there are no strict guarantees on the behavior of the diagonal elements like we get from the LLL conditions, \eqref{eq:LLLdiagonal}. Another way to think about this is to consider the geometry. The permutation based reduction only reorders a given set of basis vectors in an attempt to optimize the ordering for the search. The LLL algorithm actually finds a new, better set of basis vectors. See Figure \ref{fig:latticeBasis} for a picture of what this may look like. The LLL algorithm however has no knowledge of the input vector $y$, since we have seen that the optimal ordering for the columns of $R$ depends on $y$, it is reasonable to assume that we should be able to find better column orderings for the search process than the one which LLL gives.

The basic idea behind the proposed solution is simple, apply the LLL reduction to find a new set of basis vectors, then apply the permutation strategy given in subsection \ref{subsec:newReduction} to reorder the basis vectors. We will call this reduction strategy ``LLL+PERMU'' from here on for convenience. Unfortunately the new upper triangular matrix after the permutation reduction is applied may no longer be LLL reduced. In many cases after applying the LLL+PERMU strategy we will see a very significant decrease in search time compared to LLL reduction alone, although other times there will be no significant difference. The difference in search time between LLL and LLL+PERMU depends both on the matrices and the vector $y$ (which has a
random component). Numerical experiments indicate that when the matrix $H$ is
generated in some ways, and the standard deviation of the noise is within a
certain range, we should apply the permutation reduction, but not when $H$ is generated
in some other ways, or when the noise is too high. The performance improvement
in the search process can be quite significant on average for some practical
cases.

Denote the Babai integer point from Chapter \ref{chap:SESearch} by the vector $z_0$. This is the first point found during the SE search process. If the initial radius $\beta$ is chosen to be larger than the residual $$\left \| Rz_0 - \bar{y} \right \|_2,$$ then $\beta$ will not prune any nodes from the search tree. This means that the Babai integer point effectively defines the initial radius of the search process when the initial value chosen for $\beta$ is too large. Since the Babai integer point is usually a good estimate at the ILS solution and is cheap to compute, often we simply use an initial $\beta = \infty$ which for the search process is equivalent to using the residual from the Babai integer point. The number of nodes that the SE search visits in the search tree is strongly related to the initial radius, the ordering of the columns and the shape of the lattice (the shape of the lattice defines how many integer points are within the sphere defined by the radius). It is obvious how the initial radius and shape of the lattice relate to the search process, a smaller radius for a fixed lattice results in a smaller search space. If the lattice points are densely packed, even a small sphere radius could include many of them. To see this, again consider Figure \ref{fig:latticeBasis}. Looking at the ``Babai Point Radius'' in the figure, we can see that there are many potential solutions that lie within it, these would all be valid candidate solutions if we had started with the ``Babai Point Radius'' as our initial radius. If this initial radius had been smaller, there would be fewer potential solutions and therefore the search space would be smaller as well.

We will usually obtain different Babai integer points from the two reduction strategies, LLL and LLL+PERMU. Since we know that the search time is related to the initial radius which is often defined by the Babai integer point as mentioned above, one way to estimate whether the search process will proceed faster after permutations is to compare the initial radius in both cases. Figure \ref{fig:LLLvsPermuBabai} shows such a comparison. Looking at Figure \ref{fig:LLLvsPermuBabai}, we can see that after permutations, the Babai integer point in this case is at least as good as the one before permutation. Often it is significantly better.  This is because the permutation based reduction strategies work directly with the Babai integer point which depends on the vector $\bar{y}$. With the permutations, we are trying to maximize the probability that the Babai integer point is the ILS solution. The result is that the Babai integer point should have smaller residual in many cases.

Our goal is to complete the search process as quickly as possible. At this point we have two ILS problems which have the same optimal solution, one is LLL reduced and the other comes from applying the permutation strategy to the first, e.g., LLL+PERMU. These two problems will often have different residuals to the Babai integer point. A good strategy to try and optimize the search process is to compare the residual to the Babai integer point from these two equivalent ILS problems and perform the search process on the one with the smaller residual. Numerical results comparing this strategy to performing the search on the LLL reduced matrix each time can be found in section \ref{sec:ILSReductionResults}.

Unfortunately, there are certain special types of ill conditioned problems where the above strategy does give an improvement in search speed. Even though the new ILS problem after LLL+PERMU gives a better Babai integer point, the search process may visit a slightly larger number of nodes. Some characteristics and examples of such problems are given in section \ref{sec:ILSReductionResults}. In these cases, we may use the permutation strategy just to help set an initial search radius $\beta$ as was discussed in Chapter \ref{chap:SESearch}. Then we may perform the search on the LLL reduced matrix $R$.

We may use the residual from the Babai integer point obtained from the LLL+PERMU reduction directly as the initial $\beta$ in the LLL reduced problem. Using this residual in the search process on the LLL reduced problem when it is smaller than the residual to the Babai integer point given by LLL reduction alone is guaranteed to result in a shorter search time. From here on we will refer to this strategy as LLL+BABAI. Since the permutation reduction process is very fast compared to the search, usually the overall run time, which is the time for the search plus the time for reduction, will be faster as well. There are some cases where we obtain better results by using the matrix from the LLL+PERMU strategy in the search, and some where it is better to use the LLL+BABAI strategy. Some comparisons can be found in section \ref{sec:ILSReductionResults}.

\vsp \subsection{Numerical Experiments} \label{sec:ILSReductionResults}
In this subsection, some numerical results will be presented to compare the speed of the search process after applying the reduction strategies presented in section \ref{sec:oilsreduction}. We will compare the LLL reduction, LLL reduction plus permutations (LLL+PERMU), and the LLL reduction using the initial radius $\beta$ given by the permutation strategy (LLL+BABAI). For all of the search time results the time to find the Babai integer point is not included.

Figure \ref{fig:LLLvsPermuBabai} which was discussed in the previous subsection
shows the residual to the Babai integer point before and after the PERMU strategy has
been applied. This figure shows numerical results for $200$ matrices $H \in\mathbb{R}^{50
\times 50}$ generated randomly with their elements picked from a normal
distribution, $N(0,1)$. The elements in the vectors $x \in \mathbb{Z}^n$ were picked
uniformly on the interval $[-10,10]$ and the vector $y$ was generated as $y = Hx
+ v$, where $v \sim N(0,(0.5)^2I)$ in this case.

\begin{figure}
\centering
\includegraphics[scale=0.36]{lllvslllpermubabai.png}
\caption{LLL Reduction vs LLL+PERMU. Residual from the Babai points.}
\label{fig:LLLvsPermuBabai}
\end{figure}

Figure \ref{fig:LLLvsPermuvsBabai} shows the average search time when the three reduction strategies (LLL, LLL+PERMU, LLL+BABAI) are used to reduce some normally distributed, random square matrices $H \in \mathbb{R}^{n \times n}$ with $ n = 35:5:50$. The horizontal axis displays varying values for $\sigma$ (``Sigma'' in the figure) which is the standard deviation of the noise vector $v$. The vertical axis is the average time taken for the search process over $200$ runs where on each run a new matrix $H$, vector $x$ and noise vector $v$ were generated.

Looking at Figure \ref{fig:LLLvsPermuvsBabai}, we see that when the noise is
small, the search process is too fast for any extra reduction to make a
difference. As the noise becomes larger, LLL+PERMU and LLL+BABAI start to offer a significant
advantage in terms of search time. LLL+PERMU and LLL+BABAI perform similarly at
moderate levels of noise in this case. As the noise rises further, the search time after LLL+PERMU seems to become about the same or slightly worse than after LLL, however LLL+BABAI maintains a significant
advantage in all cases.

\begin{figure}
\centering
\includegraphics[scale=0.36]{lllvspermuvsbabai.png}
\caption{LLL Reduction vs LLL+PERMU vs LLL+BABAI. Average times over 200 runs for various problem sizes and noise levels.}
\label{fig:LLLvsPermuvsBabai}
\end{figure}

To better understand the behavior of the search process with the three reduction strategies, Figure \ref{fig:spikes} shows the search time for 200 runs on a few select problem sizes and levels of noise. Notice that when the noise is at a lower level like $0.5$, the LLL reduction has many spikes in its runtime, where LLL+PERMU and LLL+BABAI encounter these spikes much less frequently. Most of these spikes can be explained by looking at the residual to the Babai integer point. When there is a spike in the LLL time it is usually the case that the Babai integer point after LLL reduction resulted in a large residual, while the LLL+PERMU and LLL+BABAI reductions gave a smaller residual to the Babai integer point. The reason for this was explained in subsection \ref{sec:oilsreduction}. When the noise becomes larger like around $0.8$ we see large variance in the runtime from all three reduction methods. This is because when the noise becomes too large, the difference between the residuals to the Babai points between LLL and LLL+PERMU becomes smaller.

\begin{figure}
\centering
\includegraphics[scale=0.36]{spikefig.png}
\caption{LLL Reduction vs LLL+PERMU vs LLL+BABAI. Search time for 200 runs with some different sizes and noise levels.}
\label{fig:spikes}
\end{figure}

Table \ref{tab:successRates} shows how many times out of $200$ runs that the
Babai integer point was equal to the true solution $x$ for both LLL reduction and
LLL+PERMU. The label ``New'' refers to the LLL+PERMU algorithm. The vertical axis of the table shows the various levels of noise,
$\sigma$, and the horizontal gives the different problem sizes. If we take this
as a percentage it is known as the ``success rate'', it gives the probability that the Babai integer point is equal to the true parameter vector. This is an important number for GPS applications where we are required to solve OILS problems. In GPS the success rate can be estimated based on the problem. If the estimated success rate is high enough, the expensive OILS problem is not solved, the Babai integer point is simply used as an estimate for the parameter vector instead. In communications applications the success rate is generally not considered, instead there are other measure like the ``bit error rate'' and ``symbol error rate''. The success rate is shown here because we are considering OILS problems and communications applications deal with BILS problems. Here we see that the LLL+PERMU strategy can offer a significantly better success rate in all cases tested; even in cases where the search process may be slower on average if we use the permuted $R$.

\begin{table}
\caption{Success Rate (out of 200) for LLL and LLL+PERMU on various problem
sizes and levels of noise.}
\label{tab:successRates}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
$\sigma$&               \multicolumn{2}{|l|}{$35 \times 35$} & \multicolumn{2}{|l|}{$40
\times 40$} &\multicolumn{2}{|l|}{$45 \times 45$} & \multicolumn{2}{|l|}{$50
\times 50$} \\ \hline
&               LLL & New & LLL & New & LLL & New & LLL &
New \\ \hline
$0.05$&         $200$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$ 
\\ \hline        
$0.1$&           $200$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$
\\ \hline
$0.3$&           $199$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$ & $200$
\\ \hline
$0.4$&           $192$ & $200$ & $196$ & $200$ & $193$ & $200$ & $196$ & $199$
\\ \hline
$0.5$&           $178$ & $198$ & $175$ & $198$ & $175$ & $194$ & $167$ & $196$
\\ \hline
$0.7$&            $89$ & $132$ &  $75$ & $120$ & $79$  & $141$ & $89$  & $137$
\\ \hline
$0.8$&            $44$ &  $76$ &  $45$ &  $77$ & $42$  & $86$  & $44$  & $94$
\\ 
\hline
\end{tabular}
\end{table}

Looking at these numerical results, we see that when the noise is both too low, or too high there seems to be no big difference between the reduction strategies. When the noise is too low this is because the problem becomes too easy, all three reductions result in a very fast search speed. When the noise is too high the explanation is a little bit more complicated. We have performed numerical experiments which concluded that as the noise becomes high, the new upper triangular matrix given by LLL+PERMU becomes farther from a LLL reduced matrix, e.g., the Lovasz condition \eqref{eq:LLLConditions2} is less likely to be satisfied. Since we know that the properties of LLL are desirable in the search process, it makes sense that if they are lost the search may slow down. This problem still needs further investigation for a more rigorous explanation.

The next set of numerical results demonstrates a case where the permutation strategy should not be used. In Figures \ref{fig:lllvspermubabaiill} and \ref{fig:lllvspermuvsbabaiill} the matrix $H$ is generated as follows.  We construct a diagonal matrix $D = diag(d_{ii})$,  where $d_{ii} = 10^{-\frac{(i-1)*4}{(n-1)}}$. Recalling that the condition number of a matrix is the first singular value divided by the last one; a  matrix with condition number $10^4$ is constructed by forming the product of the SVD, $H = UDV^T$. Here $U$ and $V$ are the $Q$ factors from the QR decomposition of two random matrices. We generate $x$ and $y$ in the same way as the previous numerical experiments.

\begin{figure}
\centering
\includegraphics[scale=0.36]{lllvspermubabaiill.png}
\caption{LLL Reduction vs LLL + PERMU. Residual from the Babai points on ill conditioned problems.}
\label{fig:lllvspermubabaiill}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.36]{lllvspermuvsbabaiill.png}
\caption{LLL Reduction vs LLL +PERMU vs LLL + BABAI, average search time on ill conditioned problems.}
\label{fig:lllvspermuvsbabaiill}
\end{figure}

Figure \ref{fig:lllvspermubabaiill} shows the residual for the Babai integer point both before and after permutation where $H \in \mathbb{R}^{40 \times 40}$ and $\sigma = 0.5$. Comparing to Figure \ref{fig:LLLvsPermuvsBabai}, we notice that it is different. First, the residual both before and after permutation is much smaller in magnitude. The residual after permutation is still slightly better than before on average, but the difference is much less significant than Figure \ref{fig:LLLvsPermuvsBabai}. The residual to the Babai integer point after LLL reduction alone has less variance relative to LLL+PERMU compared to Figure \ref{fig:LLLvsPermuvsBabai}. Out of the $200$ runs in Figure \ref{fig:lllvspermubabaiill}, not once was the Babai integer point equal to the true solution; compared to the corresponding example in Figure \ref{fig:LLLvsPermuBabai} where $167/200$ times the Babai integer point was the correct solution before permutations were performed. 

Figure \ref{fig:lllvspermuvsbabaiill} shows the average search times for these ill conditioned problems. We see that LLL+PERMU is consistently slightly worse than the other reduction strategies. In such ill conditioned cases, using LLL+BABAI is the best strategy.

\chapter{Alternate Search Strategies} \label{chap:Searches}

In this chapter some more search strategies will be reviewed. The best first search which has been mentioned previously will then be presented in detail. In the literature there have been a few attempts to control the memory usage of the BFS, one of these will be looked at and a new one will be presented along with some numerical results.

Chapter \ref{chap:SESearch} described the most common search strategy used to solve both the OILS and BILS problems. Also, the search process was shown to be equivalent to a tree search problem where we must find the minimum cost leaf in a tree of depth $n$ with potentially exponential width.

While the depth first search corresponding to the SE algorithm is usually very fast, it is not optimal in terms of the number of nodes visited during the search process. Let $x_{i:n}^p$ denote the node in the search tree that results from fixing $x_{i:n}$ to some particular set of values, where the superscript $p$ stands for partial. We may call this a partial solution. The squared partial residual given by this partial solution is just its cost in the search tree and can be defined as $$\left \| R_{i:n,i:n}x_{i:n}^p - \bar{y}_{i:n}\right \|_2^2.$$ 

Consider the search problem only as a tree search, assume we have the nodes and edges with their weights but have no other knowledge of the problem. Then we can define an optimal algorithm for an ILS problem as one which visits only nodes in the search tree with partial residuals that are less than or equal to the optimal ILS solutions residual. Therefore all nodes visited by an optimal search process satisfy the following inequality:
\begin{equation}\label{eq:optimalSearchCond}
\left \| \bar{y}_{i:n} - R_{i:n,i:n}x_{i:n}^p \right \|_2^2 \le \left \| \bar{y} - Rx_{\scriptscriptstyle{\mathrm{ILS}}} \right \|_2^2.
\end{equation}
Here  $R$ and $\bar{y}$ come from reducing an ILS problem defined by $H$ and $y$, and $x_{\scriptscriptstyle{\mathrm{ILS}}}$ denotes the optimal ILS solution, which depending on the noise may not be the same as the original integer vector. 

Assuming no other knowledge of the search problem beyond the tree, any search algorithm must at least visit this set of nodes to guarantee there is no other leaf in the tree with a smaller residual than the one found (recall from Chapter \ref{chap:SESearch} leaves all lie at depth $n$). This means all optimal search algorithms as defined above visit the same set of nodes but possibly in a different order. The SE search presented in Chapter \ref{chap:SESearch} is not optimal since some of the nodes which it visits may not satisfy \eqref{eq:optimalSearchCond}. 

There are in fact ways to reduce the number of nodes visited so that we may visit fewer nodes than the optimal algorithms as defined above. One technique involves using lower bounds on the optimal solutions residual to prune nodes from the search tree, see e.g., \cite{StoVH08}. Unfortunately our numerical tests indicate that existing methods which use this technique incur too much per-node overhead to provide faster solutions to {\em reduced} OILS and BILS problems.

One algorithm that satisfies the requirement in \eqref{eq:optimalSearchCond} is called the ``Best First Search'', which will be referred to as BFS from here on for convienience. This algorithm has been proposed in the literature a number of times, \cite{FukMU04} and \cite{XuWZW04} are two examples, although they appear different on the surface. Later in this chapter a detailed description of the BFS will be given. For now it is enough to note that the BFS pays a price to achieve the goal of visiting the minimum number of nodes, it must permanently store each node visited during the search in memory and always keep track of the node with minimum cost (where cost is defined as in Chapter \ref{chap:SESearch}). Also, when the nodes are later accessed from memory, since they are not stored in a simple manner like the SE search, it is likely that we will incur cache penalties when accessing them. For these reasons it is not always practical to use the BFS strategy since hardware applications may have very limited memory and slow processors.

Unfortunately, some of the literature such as \cite{StuBF07} does not fully consider the overhead of finding the minimum cost node when comparing search algorithms, instead only considering limiting the memory usage of BFS while still visiting as few nodes as possible. Even though such algorithms may visit fewer nodes than an SE search, the extra time spent at each node may cause them to still run slower. In section \ref{sec:controlBFS} the algorithm from \cite{StuBF07} will be described. Comparing the number of nodes visited by two different search processes may not always be meaningful if one process spends much more time on each node than the other.

The following sections will explain the BFS, quickly overview a couple of attempts that have been made to control the memory usage of the BFS, and then propose a new idea for combining the BFS with the SE search in order to limit memory usage and decrease computational complexity.

\vsp \section{Best First Search} \label{sec:BFS}

Recall from Chapter \ref{chap:SESearch} that we will explore the search tree which has depth $n$ and each edge has a weight. Define a  nodes cost as the sum of the  weights of all edges traversed to reach the node, that is the length of the path between the node and the root. This cost is equal to the partial residual squared for fixing $x_{i:n}$ to a set of values, i.e., $$\left \| R_{i:n,i:n}x_{i:n}^p - \bar{y}_{i:n}\right \|_2^2.$$ We would like to find the leaf with minimal cost. For the BFS, we define a node's ``next best child" as the child of that node with the lowest cost that has not yet been visited. This is simple to compute; the first ``next best child'' for a node is just $\lfloor c_k \rceil$ as defined in Chapter \ref{chap:SESearch}. The second ``next best child'' for the node is $\lfloor c_k \rceil \pm 1$, such that we get the second nearest integer to $c_k$, we proceed in this way, always taking the next nearest integer to define the ``next best child''. This is the exact same strategy that the SE search uses.

The core data structure used to efficiently implement a best first search is called a priority queue. The priority queue is an abstract data structure whose basic operations are $insert(element,cost)$, and $findmin()$. The $insert(element,cost)$ operation adds an element to the queue with some real value cost. The $findmin()$ operation finds and removes the element with minimum cost from the queue. The implementation details are not particularly relevant to this thesis, but it is important to note that there are various implementations used in practice. The usual cost for the $insert(element,cost)$ operation is $\theta(1)$ and for the $findmin()$ operation $\theta(log(N))$ where $N$ is the number of elements currently in the queue. For more information on priority queues and their implementation see \cite{AhoHU83}.

As mentioned above, for an ILS application, we can quickly find the next best child of a given node, it is the node which corresponds to the next choice of $x_k$ chosen as in the SE algorithm. The ``first best child'' of a node at level $k-1$, assuming $x_{k-1:n}$ are fixed, is always the node corresponding to $x_k = \lfloor c_k \rceil$.

In the first step of the BFS, we initialize an empty priority queue and define the root node as having a cost of $0$ and being at level $n+1$. The ``next best child" of the root node corresponds to $x_n = \lfloor \bar{y}_n/R_{n,n} \rceil$ and the cost to visit this child is equal to the partial residual squared $(R_{n,n}x_n - \bar{y}_n)^2$, we call this cost the ``next best child cost". The root node will store its current ``next best child" and ``next best child cost" and then be inserted into the priority queue. Elements in this priority queue are always sorted by their ``next best child cost".

In the next step, we visit the first child of the root, $x_{n}$. First, we perform the $findmin()$ operation on the priority queue, since currently the root is the only element, it will be returned. We would like to visit the first child of the current node (which is now the root). To visit a node involves calculating the ``next best child" and the ``next best child cost", therefore we must compute these quantities for the node corresponding to $x_{n}$ (which is the first child of the root). The ``next best child" of $x_n$ will be given by $x_{n-1} = \lfloor c_{n-1} \rceil$ and its cost is the partial residual squared: $$\left \| R_{n-1:n,n-1:n}x_{n-1:n}^p - \bar{y}_{n-1:n}\right \|_2^2.$$ We then insert $x_n$ into the priority queue with its ``next best child'' and ``next best child cost''. Notice that if we expand the cost as follows, $(R_{n,n}x_n - \bar{y}_n)^2 + (R_{n-1,n}x_n + R_{n-1,n-1}x_{n-1} - \bar{y}_{n-1})^2$ then the first term in the cost is just the cost of the parent node; this means to calculate the cost for a child at level $i$, we only need to add the $(n-i+1)^{th}$ term from the sum defined by $$\left \| R_{i,i:n}x_{i:n}^p - \bar{y}_{i}\right \|_2^2$$ to the cost of the parent.

Since we are now visiting the first child of the root, we must generate the new ``next best child" for the root node, $x_n = \lfloor c_n \rceil \pm 1$ and the cost for this child (also the partial residual), $(y_n - R_{n,n}x_n)^2$. We may now insert the root back into the priority queue with the newly calculated ``next best child cost" and ``next best child''.
 
At this step, the next node to be visited will be the one at the top of the priority queue with the smallest ``next best child cost". Currently there are two nodes in the priority queue, one at level $n+1$ (the root) and one at level $n$. If the former has a smaller ``next best child cost'', we will visit the second best child of $x_{n+1}$, the new $x_{n}$ next. Otherwise we will visit the first child of the node $x_{n}$, which is at level $n-1$.

By proceeding in this way we always visit the nodes in the order of increasing partial residual. The next node we visit is always the one with the next smallest partial residual (in the whole tree) from the previous one, even if those 2 nodes are at different levels. This should be obvious, consider the case where there exists a node in the tree that we have not yet visited which has a smaller partial residual than the current node at the top of the priority queue. This means that the parent of the node in question also had a smaller partial residual than the current node and its grandparent etc. up to the root. If this were the case, however, we would be visiting that node instead of the current node, so it is a contradiction. In this way we can guarantee that the first time we find a leaf in the tree, it must be the leaf with the smallest squared residual and is therefore the solution. Also, at the point where we find a leaf, we know we have visited only and all of those nodes that have a partial residual within the radius defined by the optimal solutions residual.

It should also be noted that trivial modifications to this BFS search algorithm can be made so that it will solve the BILS problem. We just need to make sure that at every step, the value chosen for $x_k$ is within the proper bounds defined in equation \eqref{eq:boxCon}.

\vsp \section{Controlling BFS Memory Usage} \label{sec:controlBFS}

In section \ref{sec:BFS}, notice that in each step a new node is added to the priority queue, but nodes are never removed (the $findmin()$ operation removes a node, but we just update the cost and insert it back in). This means that each node visited during the search is kept permanently in the priority queue. Since the number of nodes visited can potentially be exponential in $n$ this can become a problem for two reasons. The first and most obvious reason is memory usage. The second which is often overlooked is that, as the priority queue grows, the computational complexity to maintain it is increased at each step; the computational cost for each operation is $\theta(log(N))$, where $N$ is the number of nodes visited so far. Also consider that the computational cost to visit a node in an efficient implementation of the SE algorithm is only about $2k$ flops, where $k$ is the level of the node in the tree. The BFS must still perform the operations the SE performs in addition to maintaining the priority queue.

With the rough analysis of the computational cost in the previous paragraph, it should be obvious that just because the BFS visits fewer nodes than SE, it does not mean it will necessarily be faster. In both \cite{DaiY08} and \cite{StuBF07} the authors try to find a balance between the number of nodes that we keep in memory at any given time, and visiting the smallest number of nodes possible; this section will briefly describe the approach taken in \cite{StuBF07}. Note that the algorithm described here is slightly different in that an extra parameter has been eliminated, according to the authors results, this parameter seems to always make the performance worse anyways.

Their algorithm is only slightly different from the BFS. The authors do not use the concept of a priority queue, and instead use ordered lists to find the node with minimal cost at each step. Suppose we have a list $S$ in which nodes are stored, also this list is sorted by each nodes level in the search tree (lower levels are toward the back). This list, like the priority queue in the BFS starts with only the root node. At each step, a new node is added to the back of the list by finding the node with minimal cost in $S$ (this takes $|S|$ operations) and visiting its ``next best child''. Now suppose we allow the user to specify some integer parameter $\alpha$. Instead of scanning the whole list $S$ to find the node with the lowest ``next best child cost'' to visit next, we simply look at the last $\alpha$ elements in $S$ (which correspond to the $\alpha$ lowest nodes in the tree) and make our decision based on these. Such a strategy forces the BFS to proceed down the tree much faster than it would otherwise.

Unfortunately, when a leaf is reached, we no longer have the guarantee that it is the optimal solution. We do know however that we no longer have to consider any of the last $\alpha$ nodes in $S$ since they all must have a cost greater than the cost of the leaf we had just found. We may remove the last $\alpha$ nodes in $S$, update the search radius to be the residual given by this leaf, and continue the search process. Any time a node is visited with a ``next best child cost'' that is greater than the current search radius, we may discard the last $\alpha$ nodes in $S$. When $S$ is empty we may terminate with the optimal solution.

The behavior of this algorithm is such that when $\alpha = 1$ it degrades to a SE search. When $\alpha = \infty$ it becomes a BFS. Unfortunately implementing a BFS using a list is very inefficient as will be discussed later. The authors also give some good bounds on the amount of memory this algorithm requires based on how the parameter $\alpha$ is chosen. They present some results, but do not give flop counts or CPU time, instead focusing on memory usage and the number of nodes visited.

There are a few drawbacks to this algorithm. One is that it is not clear how to implement this in practice. Consider setting the parameter $\alpha$ to a relatively high number. When $\alpha$ is higher, we will visit fewer nodes since it will be closer to the BFS, and we will discard more nodes each time we have the opportunity to discard. Unfortunately as $\alpha$ gets larger, a naive list based implementation becomes impractical. Scanning through $\alpha$ elements in a list at each step could impose significant overhead. Also consider that $S$ must be sorted by the levels of the nodes in the search tree. If we wish to add a new node which has a level larger than the smallest leveled nodes currently in $S$, we must move all of the lower nodes in the list one place to the right in memory, again this could take $\alpha$ operations. This suggests that we may want to use a priority queue based implementation for large $\alpha$ and a list based implementation for small $\alpha$. A priority queue implementation of this algorithm is not very straight forward however since we have the extra condition that the node at the top of the priority queue should have the lowest cost and a level in the tree corresponding to the parameter $\alpha$. 

With these problems, it is worthwhile to explore some other options to limit the memory usage of the BFS, while trying to decrease the amount of time spent processing each node in the tree.

\vsp \section{Combining BFS and SE Search}

Here a new method for controlling the memory usage of the BFS will be presented. Like the method mentioned in the previous section \cite{StuBF07} it is a modification to the BFS which relies on a parameter supplied by the user.

One of the main problems with the SE search occurs when elements high up in the tree (e.g. $x_n, x_{n-1} , \dots$) are initially chosen incorrectly. We must explore the entire subtree looking for the ILS solution, visiting many nodes that may have partial residuals higher than the residual given by the ILS solution. Eventually we make our way back up the tree to correct the mistake. The BFS does not have this problem since it does not restrict itself to visiting only children of the current node.

In order to avoid making expensive mistakes at high levels, we may consider dividing the search tree into two parts, cutting it at some level $\alpha$ which is to be supplied by the user. Nodes with level greater or equal to $\alpha$ will belong to one part and less than $\alpha$ another. For example, consider figure \ref{fig:treeSearch}. If we set $\alpha = 2$, the nodes labeled $1,6,2$ and $4$ would belong to the top group, while nodes $3$ and $5$ would belong to the bottom. We begin with a BFS on the top part of the tree, when the BFS reaches a node at level $k=\alpha$, it uses an SE search to find the optimal solution in the subtree of that node. Suppose we are at some node in the BFS which is given by fixing $x_{k:n}$ to some particular values (so we are at level $k=\alpha$), then we solve the subproblem  given by the following using the SE search.

Let $$\hat{y} = \bar{y}_{1:k-1} - R_{1:k-1,k:n}x_{k:n},$$ where $\bar{y}$ and $R$ come from one of the reduction strategies described in Chapter \ref{chap:Reduction}. Note that the search process in Chapter \ref{chap:SESearch} also combines the known part of $x$ with $R$ and $\bar{y}$ as it progresses in the same way. Then assuming that $x_{k:n}$ are fixed, we would like to find the $x_{1:k-1}$ to minimize the following: $$\min_{x_{1:k-1} \in {\mathbb{Z}^{k-1} }}  \| \hat{y} - R_{1:k-1,1:k-1}x_{1:k-1} \|_2.$$

Call the estimate $x_{1:n}$ which is given by combining the solution to the above subproblem ($x_{1:k-1}$) with the nodes found by the BFS ($x_{k:n}$) the ``current best integer point''. We can continue the BFS until we reach another node at level $\alpha$. If at some point during the BFS the next node in the priority queue has a ``next best child cost'' greater than the cost of the current best integer point found by a previous SE search, we may terminate the search process knowing that the current best integer point is the optimal ILS solution. When the BFS eventually visits another node at level $\alpha$, we again perform an SE search to find the optimal solution in the subtree of this node, however this time we may use a sphere constraint to terminate the SE search early i.e., we can set a good initial $\beta$. We can use the current best integer point to set the initial $\beta$ for the search on this new subproblem. To do this we simply set the initial $\beta$ for this new search equal to the cost of the current best integer point minus the cost of the node at level $\alpha$ we are visiting, where we will start our SE search. This ensures that any leaf found by the current SE search must have a lower cost than the current best integer point. If this SE search finds a leaf, we may update the current best integer point to correspond to this new leaf. It is possible that the current SE search will terminate without finding a leaf, this means that the corresponding subproblem does not yield a better solution than the current best integer point. Continuing this process, eventually we will reach a point where the next node to be visited in the BFS has a cost greater than the current best integer point, at this point we may terminate with the optimal solution.

With this approach and an appropriate choice of the parameter $\alpha$, we can hope to avoid making expensive mistakes at higher levels in the search tree. We also will incur a much smaller overhead to maintain the priority queue, since it will not contain nearly as many nodes if $\alpha$ is chosen to be near $n$, and with the smaller priority queue, there will be less memory usage. 

There are some useful tricks which can be used in conjunction with the combined search process to improve performance. Consider that the first SE search will very often be the one that costs the most in terms of run time. This is because we do not have a good value for an initial $\beta$ to use as a stopping condition and must solve the subproblem completely, even if the subproblem may be very ill conditioned. The best case occurs when the first subproblem we solve yields the ILS solution, in this case we get the smallest possible sphere radius to use in subsequent SE searches as a stopping criteria, and therefore we can prune more nodes from the search tree. To make this more likely to occur, we may delay solving the first subproblem as follows. When we visit a node at level $\alpha$, we can decide not to solve the subproblem stemming from this node right away. Instead we put the node into a list of ``unsolved subproblems'' and continue the BFS. One example of when we may not want to solve a subproblem immediately is when the Babai integer point given by that subproblem has a large residual. Eventually we need to perform the SE search to find a leaf. Once a suitable node is visited by the BFS (one with a relatively low residual to the Babai integer point is a good choice), we perform the SE search and find a leaf in the tree corresponding to a potential solution. Now we simply go back to the list of ``unsolved subproblems'' and solve them all, using the residual from the previous solution as an initial sphere radius. The unsolved problems that had Babai points yielding high residuals are now likely to terminate much more quickly because not many nodes should be within the current sphere radius if the subproblem that we chose to solve gave a good solution. This strategy has not been explored in full detail and still requires some future research on the best way to utilize it.

Finally, it should be mentioned that while this combined BFS and SE search strategy has been described in a way to solve OILS problems, it is also applicable to BILS problems. Only minor modifications must be made in order to ensure the elements of $x$ satisfy the constraints defined in \eqref{eq:boxCon}.

The advantage of this new algorithm over the one presented in \cite{StuBF07} which was described in section \ref{sec:controlBFS} is that the algorithm in \cite{StuBF07} has no clear and efficient implementation. The implementation you should use depends how you would like to choose the parameter $\alpha$. Even for a fixed $\alpha$ it is not always obvious exactly how the algorithm should be implemented efficiently. The new algorithm presented here does not have this problem. The implementation is straight forward and should always use a priority queue, no matter how the parameter is to be chosen.

Some strategies for choosing the parameter $\alpha$ based on the input problem were tested. Unfortunately choosing $\alpha$ based on properties of the problem has proven to be somewhat difficult, although it can be chosen reliably through numerical experiments as will be shown in the next section. How to choose $\alpha$ for a given problem will be considered in future research.

\vsp \section{Numerical Testing Results} \label{sub:SearchResults}
In this section, numerical results will be given to compare the SE search, BFS and the combination of the two. To this point, there has not been much said about how to choose the parameter $\alpha$. We need to choose $\alpha$ low enough to avoid mistakes at high levels in the tree, but high enough so that we get the advantages of the SE search being faster at processing each node. Note that when $\alpha = n$, we have a pure SE search, and when $\alpha = 1$, the combined search strategy degrades to a BFS. 

For the first set of numerical results in this section, the matrices are simply generated randomly with elements drawn from a normal distribution, $N(0,1)$. The integer vectors $x$ have elements uniformly distributed in the interval $[-10,10]$ and the vectors $y$ are generated as $y = Hx+v$ where $v$ is normally distributed with some standard deviation $\sigma$. All matrices in this section will be LLL reduced before the search process begins.

Figure \ref{fig:searchsurface} shows runtime results as the parameter $\alpha$ varies for $50$ different randomly generated OILS problems of size $50 \times 50$. Here the noise vector has $\sigma = 0.8$, which is quite high. The noise was chosen to be high so that the runtime would be significant, when the noise is too low all of the search algorithms will be very fast. The x-axis which varies from $1 , \dots, 50$ represents the choices of the parameter $\alpha$. The y-axis gives the search time in seconds, and each line down the z-axis represents one of the $50$ random problems. On the far right of the graph where $\alpha$ is small, we have a search process equivalent to the BFS and toward the left as $\alpha$ gets larger it is equivalent to the  SE search. Notice that on both ends there are sometimes spikes in the run time, while if we follow the line defined by $\alpha = 40$, we see that there is a consistent trough in the runtime. This is interesting and tells us that for this particular type of problem, it is best to choose the parameter $\alpha$ somewhere around $40$. Also notice that choosing the parameter slightly incorrectly is usually not a problem, there is a fairly wide range of good choices. 

Figure \ref{fig:searchAverage} shows the mean search time over all 50 runs for each setting of the parameter $\alpha$. Notice that the minimum occurs around $\alpha = 40$. This confirms our choice in the previous paragraph.

We know the search behaves fairly consistently as the parameter $\alpha$ changes, and there is a wide range for good choices of the parameter; to choose $\alpha$ for some new type of problem generated in some other way, it is recommended to simply run it a few times with various parameter settings and choose the one that performs best.

\begin{figure}
\centering
\includegraphics[width=6in, height=4in]{searchsurface.png}
\caption{Run time results for the combined search process with varying parameter $\alpha$ on random problems.}
\label{fig:searchsurface}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.36]{averageSearchTimes.png}
\caption{Average run time results for the combined search process with varying parameter $\alpha$ on random problems.}
\label{fig:searchAverage}
\end{figure}

Figure \ref{fig:illConditionedSearch} shows some numerical results from running the new combined search process on ill conditioned problems. These problems were generated in the same way as the ill-conditioned problems in section \ref{sec:ILSReductionResults} and Figures \ref{fig:lllvspermubabaiill} and \ref{fig:lllvspermuvsbabaiill}. The problem size here is $35 \times 35$ and the noise vector was generated with $\sigma = 0.3$. Here we choose the noise vector to have a smaller standard deviation because these ill conditioned problems take much longer to solve. The best results are obtained when $\alpha$ is near $n$ which means we start the SE search process immediately. These examples illustrate that the best first search is not always better than the SE search even though it visits fewer nodes in the search tree.

Some reasons for this can be found by examining the number of nodes visited during the search process. First, we look at the ratio $\frac{\textit{Nodes visited by BFS}}{\textit{Nodes visited by SE Search}}$. For the previous set of numerical results where the new combined search achieved significantly better results, this ratio was $0.44$ on average, for the new set of ill-conditioned problems, it is $0.8$ on average. This means that there is not as big a difference between the best first search and SE search as far as the number of nodes visited is concerned. The fact that the BFS explores $20\%$ fewer nodes on average is not enough to offset the extra time spent on each node.

Also with these ill conditioned problems, the run times are all on a similar scale. If we look at the variances of the run times between the two figures, the variance for the ill-conditioned problems is only about $170$, while the randomly generated problems have a runtime variance of nearly $5000$. These variances were calculated as the variance of all of the run times in each figure. If we look at Figure \ref{fig:searchsurface} we can see that most of the time the SE search is fast, it is only $8$ or $9$ out of the $50$ randomly generated problems where it becomes very slow. This small number of bad cases makes it much worse on average.
 
\begin{figure}
\centering
\includegraphics[width=6in,height=4in]{illConditionedSearch.png}
\caption{Run time results for the combined search process on ill conditioned problems with varying parameter $\alpha$.}
\label{fig:illConditionedSearch}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.36]{averageSearchTimesIll.png}
\caption{Average run time results for the combined search process with varying parameter $\alpha$ on ill conditioned problems.}
\label{fig:searchAverageIll}
\end{figure}

\chapter{Conclusions and Future Work} \label{chap:Conclusion}

This thesis focused on improving both the reduction and search processes in solving OILS and BILS problems. In Chapter \ref{chap:Reduction} it was shown that two effective reduction strategies in the literature for the BILS problem are theoretically equivalent. Using the knowledge of their equivalence a new reduction strategy was proposed that uses ideas from each to compute the same solution more efficiently and in a numerically stable manner. Also in Chapter \ref{chap:Reduction}, it was shown that under some practical conditions, this new reduction process can be used to improve the performance of the OILS search process as well. This improvement seems to be a result of the Babai integer point given by the new reduction having both a smaller residual and higher success rate (defined as the probability that the Babai integer point is equal to the true solution $x$). Since the Babai integer point has a tendency to have smaller residual and higher success rate after applying these permutations, we can conclude that for both OILS and BILS problems, applications which use the Babai integer point as an estimate to the ILS solution should apply this permutation strategy in order to improve their success rate.

Some more investigation into the new reduction algorithm is still needed. The new reduction presented in subsection \ref{subsec:newOILS} has two phases, first we do the LLL reduction which uses IGTs and permutations together to achieve some properties in the upper triangular matrix $R$. The second phase applies another set of permutations to $R$ to try and achieve some different properties. We would like to investigate reduction algorithms that try to consider both sets of these properties in a single phase instead of having two disconnected phases. This will require a deeper understanding of exactly why the new LLL+PERMU reduction gives good results.

In Chapter \ref{chap:Searches} the ``Best First Search'' approach was described. The BFS has the property that it visits the fewest nodes of any search algorithm. Unfortunately the BFS is not always ideal because of its high memory usage and complexity per node visited in the search tree. Previous attempts to improve upon the BFS were mentioned and a new one introduced. 

A previous attempt at improving the BFS, \cite{StuBF07}, involved limiting the amount of memory used by the BFS using a user supplied parameter. Unfortunately it is difficult to implement efficiently which leads to a high per-node complexity. The new algorithm introduced in this thesis tries to address both the problem of memory usage and high complexity per node while remaining easy to implement. Run time results were given comparing this new algorithm to the SE search process and the improvement was shown to be significant for some practical cases. 

It would be interesting to see if some theory could be worked out for how to choose the parameter $\alpha$ in the new search method. It is currently unclear how to choose it given a particular ILS problem. Since the optimal value for the parameter seems somewhat stable for fixed problem type, size and levels of noise, this may be possible. Also it is possible that we could do away with the parameter and dynamically decide to stop doing a BFS and start the SE search when some condition is met. For the time being it is recommended to simply compare a small number of test runs with various settings of the parameter in order to choose it for some specific applications.



\bibHeading{References}
\bibliography{../ILS}
\bibliographystyle{plain}

\index[abbr]{LS@LS: Least Squares}
\index[abbr]{ILS@ILS: Integer Least Squares}
\index[abbr]{OILS@OILS: Ordinary Integer Least Squares}
\index[abbr]{BILS@BILS: Box-constrained Integer Least Squares}
\index[abbr]{SE@SE: Schnorr-Euchner enumeration algorithm}


\printindex[keylist]{Index}{Index}{}
\printindex[abbr]{KEY TO ABBREVIATIONS}{KEY TO ABBREVIATIONS}{}

\end{document}


 






